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ABSTRACT 

A self-contained  analysis  for  arbitrary  circulation  controlled 
airfoils  in  Incompressible  flow  is  developed.  The  analysis  predicts 
the  blowing  slot  conditions  required  to  produce  a specified  lift  co- 
efficient on  a given  airfoil  with  given  free  stream  conditions. 

An  iterative  procedure  is  used  to  find  the  blowing  slot  con- 
ditions that  allow  the  Thwaltes  condition  of  constant  pressure  in 
the  separated  region  to  be  satisfied.  With  the  input  given,  a pot- 
ential flow  analysis  is  performed  using  the  Theodorsen  method. 

Boundary  layer  analyses  for  the  lower  and  upper  surfaces  then  yield 
the  separation  pressure  on  the  lower  surface  and  the  boundary  layer 
properties  at  the  slot  on  the  upper  surface.  The  flow  is  initially 
laminar  and  usually  becomes  turbulent.  The  Cebeci,  Smith  finite  dif- 
ference method  is  used  and  an  eddy  viscosity  model  is  used  for  tur- 
bulent flow.  Blowing  slot  values  are  assumed  and  a turbulent  wall 
jet  analysis  is  performed  to  determine  the  wall  pressure  at  separa- 
tion on  the  upper  surface.  If  the  two  separation  pressures,  upper 
and  lower,  do  not  agree,  new  slot  values  are  assumed  and  the  wall  jet 
analysis  is  repeated.  Wall  jet  calculations  Include  curvature  effects 
on  the  jet  and  external  flows,  and  approximate  corrections  for  large 
pressure  gradients.  An  eddy  viscosity  model  is  used  but  is  modified 
to  allow  a negative  shear  stress  at  the  velocity  maximum.  A fully 
developed  laminar  or  turbulent  channel  flow  is  used  for  the  slot  pro- 
file. A finite  difference  method  based  on  the  Keller,  Cebeci  method 
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is  used  for  wall  jet  calculations. 

Calculations  are  performed  and  compared  to  Kind's  experimental 
data.  Results  are  in  good  agreement  for  the  bloving  coefficient  and 


local  flow  details 
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1.  INTRODUCTION 

Interest  in  STOL  flight  has  led  to  an  interest  in  circulation 
controlled  airfoils.  A circulation  controlled  airfoil  has  a bluff 
trailing  edge  with  a tangential  blowing  slot  located  slightly  up- 
stream of  the  trailing  edge.  Experiments  by  Rind  (Refs.  4,  5,  15) 
and  later  by  others  have  shown  that  such  an  airfoil  can  produce  re- 
latively high  lift  coefficients  for  relatively  low  blowing  rates  from 
the  slot.  Circulation  controlled  airfoils  thus  provide  a relatively 
simple  and  efficient  method  of  obtaining  STOL  flight. 

The  principle  of  a circulation  controlled  airfoil  is  based  on 
the  fact  that  there  are  any  number  of  valid  potential  flow  solutions 
(each  with  a different  value  of  the  circulation,  and  therefore  lift) 
for  an  airfoil  with  a bluff  trailing  edge.  In  a real  (viscous)  fluid, 
the  circulation  that  will  actually  develop  depends  on  the  separation 
characteristics  on  the  upper  and  lower  surfaces.  By  introducing 
blowing  on  the  upper  surface,  separation  on  the  upper  surface  is  de- 
layed. As  the  blowing  is  increased,  the  upper  surface  separation 
point  moves  further  back  around  toward  the  lower  surface.  As  this 
occurs,  the  potential  flow  rear  stagnation  point  also  moves  down,  and 
thus  the  front  stagnation  point  and  the  circulation  are  altered.  The 
blowing  therefore  causes  the  flow  over  the  entire  airfoil  to  be  alter- 
ed. A relatively  small  change  in  the  location  of  the  stagnation 
point  results  in  a large  change  in  the  circulation,  and  therefore  in 
the  lift.  Since  the  blowing  is  used  only  to  delay  separation,  the 
amount  required  is  small  compared  to  that  required  for  a jet  flap. 
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The  blowing  is  effective  in  delaying  separation  because  of  the 
Coanda  effect.  This  effect  (Ref.  2)  refers  to  the  ability  of  a plane 
wall  jet  to  follow  the  contour  of  a convex  curved  boundary  adjacent 
to  it.  Downstream  of  the  blowing  slot,  a wall  jet  develops  on  the 
convex  surface  of  the  airfoil.  Because  of  the  Coanda  effect,  the 
wall  jet  tends  to  remain  attached  to  the  surface.  Since  the  wall  jet 
is  considerably  more  energetic  than  the  conventional  boundary  layer 
that  would  exist  if  the  wall  jet  were  not  present,  the  wall  jet  can 
resist  separation  for  a significantly  longer  distance. 

The  flow  conditions  on  a circulation  controlled  airfoil  are  in- 
dicated in  Fig.  1.  Beginning  at  the  front  stagnation  point,  a bound- 
ary layer  develops  along  the  lower  surface  and  eventually  separates. 
Also  beginning  at  the  front  stagnation  point,  a boundary  layer  de- 
velops along  the  upper  surface  and  eventually  reaches  the  blowing 
slot.  Both  the  upper  and  lower  surface  boundary  layer  are  initially 
laminar  and  then  usually  become  turbulent.  A turbulent  wall  jet  de- 
velops downstream  of  the  blowing  slot  and  eventually  separates.  A 
typical  wall  jet  development  is  illustrated  in  Fig.  2.  A separation 
region  exists  between  the  wall  jet  separation  point  and  the  lower  sur- 
face boundary  layer  separation  point. 

For  a given  airfoil  with  given  free  stream  conditions,  the  upper 
surface  separation  point  is  determined  by  the  amount  of  blowing  used. 

A change  in  the  airfoil  shape  or  free  stream  conditions  also  results 
in  a change  in  the  separation  point,  so  the  circulation  is  a function 
of  the  airfoil  shape,  free  stream  conditions,  and  blowing  rate. 
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The  purpose  of  the  present  analysis  is  to  obtain  a self- 
contained  procedure  to  predict  the  performance  of  arbitrary  circul- 
ation controlled  airfoils.  An  outline  of  the  procedure  is  given  in 
the  next  section. 
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2.  OUTLINE  OF  THE  THEORY 

The  purpose  of  the  present  analysis  is  to  predict  the  blowing 
slot  conditions  required  to  produce  a specified  lift  coefficient 
for  a given  circulation  controlled  airfoil  with  given  free  stream 
conditions.  This  section  presents  an  outline  of  the  calculation 
procedure  and  general  considerations  involved.  The  following  three 
sections  present  details  of  the  calculations  along  with  a review  of 
the  pertinent  literature. 

Input  for  the  present  analysis  consists  of  the  airfoil  geometry, 
the  angle  of  attack  oC  , the  free  stream  Reynolds  number  , and 

the  prescribed  lift  coefficient  . The  airfoil  geometry  includes 
the  location  and  thickness  of  the  blowing  slot,  which  is  assumed  to 
be  tangential  to  the  airfoil  surface. 

With  the  input  given,  a potential  flow  analysis  is  performed. 

The  potential  flow  analysis  provides  the  location  of  the  front  stag- 
nation point,  the  rear  (potential  flow)  stagnation  point,  and  the 
velocity  and  pressure  distribution  around  the  airfoil.  The  velocity 
and  pressure  distributions  are  then  converted  to  the  form  required 
for  a boundary  layer  analysis. 

A boundary  layer  analysis  is  then  performed  for  the  lower  sur- 
face of  the  airfoil.  The  analysis  starts  at  the  front  stagnation 
point  and  proceeds  downstream  until  separation  occurs.  With  the  lo- 
cation of  the  separation  point  known,  the  pressure  coefficient  at  sep- 
aration on  the  lower  surface  Cp is  known  from  the  potential 

rSeP|_ 

flow  analysis.  The  boundary  layer  flow  is  initially  laminar,  but 
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usually  becomes  turbulent  before  separation  occurs.  For  low  Reynolds 
numbers  the  flow  may  remain  laminar,  but  In  either  case  the  analysis 

determines  Cv 

rsepL  • 

A boundary  layer  analysis  is  also  performed  for  the  upper  sur- 
face of  the  airfoil.  The  analysis  starts  at  the  front  stagnation 
point  and  proceeds  downstream  to  the  blowing  slot.  In  general,  the 
last  point  at  which  boundary  layer  calculations  are  performed  is 
slightly  upstream  of  the  blowing  slot;  however,  the  boundary  layer 
properties  at  this  point  are  used  for  the  boundary  layer  properties 
at  the  blowing  slot.  The  upper  surface  boundary  layer  is  also  init- 
ially laminar,  but  becomes  turbulent  before  it  reaches  the  blowing 
slot.  Because  of  the  prevailing  pressure  gradients,  the  upper  sur- 
face boundary  layer  will  normally  become  turbulent  even  for  relativ- 
ely low  Reynolds  numbers.  If  the  boundary  layer  separates  upstream 
of  the  blowing  slot,  calculations  are  terminated.  If  the  boundary 
layer  remains  attached  at  the  slot,  the  analysis  provides  the  bound- 
ary layer  properties  at  the  slot. 

A turbulent  wall  jet  analysis  is  then  performed.  The  analysis 
begins  at  the  blowing  slot  and  proceeds  downstream  until  separation 
occurs.  When  separation  occurs,  the  wall  pressure  coefficient  at 

separation  on  the  upper  surface  Cp  is  known  from  the  wall  jet 

SeP  jj 

analysis.  The  correct  blowing  slot  conditions  are  those  conditions 
that  result  in  a constant  wall  pressure  in  the  separated  region, 
i.e.,  Cp  * . This  is  known  as  the  Thwaites  condition 

(Ref.  10),  and  is  also  experimentally  justified.  Since  it  is  nec- 
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essary  to  know  the  blowing  slot  conditions  in  order  to  perform  the 
wall  jet  calculations  to  determine  Cps_p  , an  iterative  proced- 
ure is  required.  Blowing  slot  conditions  are  assumed  and  these  con- 
ditions, along  with  the  known  upstream  boundary  layer  properties, 
provide  the  starting  conditions  for  the  turbulent  wall  jet  analysis. 
Wall  jet  calculations  are  then  performed  from  the  blowing  slot  down- 
stream until  separation  occurs,  yielding  a value  of  Cp?ep  u . In 
general,  the  calculated  value  of  will  not  equal  the  pre- 
viously determined  value  of  C.p  . New  blowing  slot  conditions 

are  then  assumed  and  the  wall  jet  calculations  are  repeated  until 
£ t,  and  C. p agree  to  within  a prescribed  tolerance. 

Then  the  Thwaltes  condition  is  satisfied,  and  the  assumed  blowing 
slot  conditions  are  the  correct  ones  for  the  prescribed  input  con- 
ditions. 

The  procedure  described  above  constitutes  a self-contained 
analysis  of  circulation  controlled  airfoils,  provided  the  calcula- 
tions can  be  performed  for  arbitrary  airfoil  shapes  without  the  use 
of  experimental  data.  Well  developed  and  tested  general  methods  are 
available  for  the  potential  flow  and  boundary  layer  calculations. 

The  wall  jet  region  is  one  of  high  complexity.  Theories  for  this 
region  are  presently  semi-empirical  relying  as  they  do  on  experimen- 
tal data  such  as  the  surface  pressure  distribution.  A major  effort 
in  this  work  Involved  the  generation  of  a self-contained  analysis 
for  the  wall  jet  region. 
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The  potential  flow  is  calculated  by  the  Theodorsen  method  (Ref. 
7)  of  conformal  transformation.  This  is  an  exact  method  and  applies 
to  arbitrary  airfoil  shapes.  It  can  also  be  used  without  difficulty 
In  the  present  application  where  the  lift  coefficient  is  prescribed. 
Details  of  the  method  are  given  in  Section  3. 

Boundary  layer  calculations  are  performed  using  the  Cebeci, 

Smith  finite  difference  method  (Ref.  8).  For  turbulent  flow,  the 
Reynolds  stresses  are  evaluated  using  an  eddy  viscosity  coefficient, 
while  for  laminar  flow  the  eddy  viscosity  coefficient  is  set  equal 
to  zero.  The  method  provides  an  exact  numerical  solution  for  lam- 
inar flow,  and  has  been  found  to  be  one  of  the  better  methods  for 
turbulent  flow  (Ref.  13).  Transition  from  laminar  to  turbulent  flow 
is  assumed  to  occur  when  the  momentum  thickness  Reynolds  number  is 
640  for  a favorable  pressure  gradient  or  320  for  an  adverse  pressure 
gradient.  Caster's  (Ref.  9)  experimental  relations  for  the  bursting 
of  short  laminar  separation  bubbles  are  also  included  in  the  tran- 
sition cirteria  if  the  flow  is  approaching  laminar  separation.  If 
trip  wires  are  used , transition  is  assumed  to  occur  at  the  location 
of  the  trip  wire  (if  transition  has  not  already  occurred  upstream  of 
the  wire).  Separation  for  either  laminar  or  turbulent  flow  is  de- 
termined by  the  wall  shear  stress  becoming  zero.  Details  of  the 
method  are  given  in  Section  4. 

There  are  several  difficulties  that  complicate  the  wall  jet  an- 
alysis. The  presence  of  the  upstream  boundary  layer  means  that  there 
is  always  a relative  minimum  in  the  velocity  profile  for  some  dis- 
tance downstream  of  the  blowing  slot.  For  typical  circulation  con- 
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trolled  airfoil  applications,  however,  the  velocity  minimum  is  only 
significant  for  a relatively  short  distance  downstream  of  the  slot. 
The  wall  jet  profile  of  most  significance  for  circulation  control 
applications  is  therefore  one  with  a velocity  maximum  but  no  minimum. 
The  inner  portion  (below  the  velocity  maximum)  of  such  a profile  has 
features  of  a conventional  turbulent  boundary  layer,  while  the  outer 
portion  (above  the  velocity  maximum)  has  features  of  a free  turbulent 
jet.  These  analogies  are  not  exact,  however,  because  the  two  regions 
develop  simultaneously.  Eddy  viscosity  models  seem  to  give  reason- 
able shear  stress  results  for  most  of  the  inner  and  outer  regions, 
but  in  general  a somewhat  arbitrary  fairing  is  required  to  obtain  a 
continuous  eddy  viscosity  distribution.  This  is  because  the  eddy 
viscosity  in  the  outer  region  tends  to  be  much  larger  (a  difference 
of  a factor  of  10  is  not  uncommon) . Of  perhaps  more  fundamental 
importance  is  the  fact  that  an  eddy  viscosity  model  does  not  seem  to 
be  accurate  near  the  velocity  maximum.  Experiments  show  that  the 
shear  stress  near  the  velocity  maximum  is  usually  negative,  with  a 
magnitude  comparable  to  the  wall  shear  stress.  However,  an  eddy  vis- 
cosity approach  gives  a zero  shear  stress  at  the  velocity  maximum. 

In  the  case  of  a flow  with  curvature  an  eddy  viscosity  approach  gives 
a nonzero  shear  stress  at  the  velocity  maximum,  but  the  magnitude  of 
the  shear  stress  is  nowhere  near  large  enough  to  agree  with  experi- 
ment. 

For  circulation  control  applications,  additional  direct  and  in- 
direct difficulties  arise  because  the  streamlines  have  a small  radius 
of  curvature.  This  results  from  a small  surface  radius  of  curvature 
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plus  the  Coanda  effect.  A conventional  order  of  magnitude  analysis 
shows  that  the  radius  of  curvature  of  a streamline  is  ( 7Z.  f-  / ) » 

where  £ is  the  surface  radius  of  curvature  and  y is  the  normal 
coordinate  of  the  streamline.  The  normal  pressure  gradient  is  there- 
fore no  longer  negligible  and  the  pressure  varies  across  the  jet. 

From  the  irrotational  flow  results , which  must  also  be  corrected  for 
curvature  effects,  the  velocity  and  pressure  are  known  at  the  outer 
edge  of  the  jet.  However,  the  pressure  in  the  jet,  including  the 
wall  pressure,  must  be  calculated  during  the  wall  jet  solution. 

These  direct  curvature  effects  can  be  handled  without  difficulty 
by  retaining  appropriate  curvature  terms  in  the  equations;  however, 
the  indirect  effects  of  curvature  are  more  of  a problem.  The  most 
significant  indirect  effect  is  that  large  potential  flow  velocity 
and  pressure  gradients  are  associated  with  a small  surface  radius  of 
curvature.  As  a result  of  the  sustained  severe  adverse  pressure 
gradient,  the  wall  jet  thickens  rapidly,  and  a streamline  radius  of 
curvature  is  no  longer  equal  to  ( ^ -/-/).  If  this  effect  is  not 
adequately  accounted  for,  the  calculated  wall  pressure  will  not  be 
accurate,  so  separation  cannot  be  determined  accurately.  Retaining 
additional  terms  in  the  y momentum  equation  can  account  for  this 
effect.  Unfortunately,  however,  the  significant  additional  term  re- 
quired cannot  be  handled  exactly  within  boundary  layer  theory  since 
it  involves  second  order  streamwise  derivatives.  Over  the  region 
where  the  above  correction  is  required,  an  analogous  correction  to 
the  potential  flow  velocity  is  required  to  obtain  an  accurate  value 
for  the  outer  edge  of  the  jet. 
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Curvature  also  tends  to  Increase  the  shear  stress  in  the  outer 
portion  of  the  jet,  and  tends  to  decrease  it  in  the  inner  portion. 
Corrections  are  available  to  account  for  this,  but  additional  uncer- 
tainties are  introduced  since  the  outer  turbulence  spilling  over  in- 
to the  inner  region  tends  to  offset  the  decrease  in  the  inner  region. 
A final  difficulty,  which  is  not  inherent  but  is  common,  is  that 
many  circulation  controlled  airfoils  have  a discontinuous  radius  of 
curvature  at  the  blowing  slot.  This  results  from  truncating  the  air- 
foil at  the  blowing  slot  and  adding  a circular  trailing  edge  in 
place  of  the  truncated  portion.  In  general,  it  is  possible  to  do 
this  and  still  maintain  a continuous  shape  and  slope,  but  the  radius 
of  curvature  is  normally  discontinuous.  For  calculations,  the 
radius  can  be  faired  in,  but  it  seems  one  fairing  is  appropriate  for 
the  wall  jet  equations  while  another  is  appropriate  for  the  outer 
edge  conditions. 

The  basic  formulation  of  the  wall  jet  problem  consists  of  the 
boundary  layer  equations  for  the  case  where  the  surface  radius  of 
curvature  is  the  same  order  of  magnitude  as  the  thickness  of  the 
viscous  region.  The  turbulent  Reynolds  stresses  are  evaluated  in 
terms  of  an  eddy  viscosity  coefficient.  This  basic  formulation  is 
altered  to  allow  the  shear  stress  to  be  negative  at  the  velocity  max- 
imum, and  to  correct  the  y momentum  equation  for  streamline  curva- 
ture effects  that  are  not  adequately  described  by  conventional 
boundary  layer  approximations.  These  alterations  were  found  to  be 
necessary  in  order  to  obtain  useful  results,  although  the  need  for 
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them  could  probably  be  reduced  by  use  of  experimental  data  In  the 
calculations.  The  shear  stress  modification  applies  near  the  velo- 
city maximum  and  is  based  on  an  extension  of  mixing  length  theory. 

The  curvature  modification  uses  average  values  to  approximate  a term 
that  is  otherwise  unmanageable.  Both  of  these  modifications  are  pre- 
liminary, but  seem  to  be  workable.  To  obtain  the  outer  boundary  con- 
ditions, the  potential  flow  velocity  is  modified  by  a standard  cur- 
vature correction  along  with  an  additional  correction  similar  to 
that  used  in  the  wall  jet.  The  pressure  is  obtained  from  the  vel- 
ocity and  total  pressure,  where  the  total  pressure  allows  for  the 
remains  of  the  upstream  boundary  layer.  This  procedure  is  used  for 
the  outer  edge  conditions  downstream  of  the  point  where  the  velocity 
minimum  essentially  disappears.  Upstream  of  that  point,  the  wall  jet 
calculations  include  the  relative  minimum  velocity,  and  the  outer 
edge  conditions  apply  at  a point  corresponding  to  the  outer  edge  of 
the  upstream  boundary  layer.  Outer  edge  conditions  for  this  region 
are  taken  as  the  uncorrected  potential  flow  velocity  and  pressure. 

The  reasons  for  using  different  conditions  for  the  two  regions  are  to 
allow  for  the  curvature  discontinuity  that  typically  occurs  at  the 
slot,  and  to  switch  smoothly  from  a profile  with  a relative  minimum 
velocity  to  one  without  a minimum. 

The  wall  jet  equations  are  solved  by  a finite  difference  method 
based  on  the  Keller,  Cebecl  method  (Ref.  14).  Details  are  given  in 
Section  5.  Although  most  wall  jet  analyses  use  integral  methods,  a 
finite  difference  method  was  selected  since  it  offered  some  potential 
advantages,  such  as  being  able  to  determine  separation  by  the  wall 


shear  stress  becoming  zero,  having  the  possibility  of  more  flexabil- 
ity  in  handling  profiles  with  a minimum  velocity,  and  being  able  to 
start  at  the  slot  without  requiring  a separate  mixing  region  analysis. 
A disadvantage  is  that  more  detailed  information  is  required,  but  the 
more  elaborate  integral  methods  require  almost  the  same  information. 
Since  the  wall  jet  calculations  are  repeated  during  the  iterative 
procedure,  a relatively  quick  and  efficient  procedure  is  required. 
Also,  the  shape  of  a wall  jet  profile  can  require  two  to  three  times 
as  many  grid  points  as  a conventional  boundary  layer  in  order  to  ob- 
tain comparable  accuracy  in  the  solution.  Therefore,  a method  that 
is  capable  of  reasonable  accuracy  with  a relatively  coarse  grid  is 
helpful.  The  Keller,  Cebeci  method  satisfies  these  requirements  and 
is  relatively  easy  to  use  and  was  therefore  selected. 
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3.  POTENTIAL  FLOW  ANALYSIS 


The  potential  flow  is  found  by  the  Theodorsen  method  of  con- 
formal transformation.  This  Is  a direct  method  for  finding  the  po- 
tential flow  over  a given  arbitrary,  closed  two-dimensional  body. 

The  method  is  exact  although  in  practice  solutions  must  be  obtained 
numerically.  The  Theodorsen  method  is  reviewed  herein  in  sufficient 
detail  so  that  the  computer  program  in  the  Appendix  can  be  under- 
stood. Additional  details  may  be  found  in  Ref.  7. 

It  is  well  known  from  potential  flow  theory  that  the  flow  over 
a given  airfoil  can  be  found  by  conformally  transforming  the  airfoil 
shape  into  a circle.  The  Theodorsen  method  transforms  the  given  air- 
foil shape  into  a circle  with  the  aid  of  an  Intermediate  transforma- 
tion as  shown  in  Fig.  3.  The  reason  for  using  the  intermediate  tran- 
sformation and  for  taking  the  airfoil  coordinate  system  as  shown  is 
that  convergence  of  the  final  transformation  is  enhanced. 

The  first  transformation  maps  the  airfoil  shape  ( 2 plane) 

into  a "nearly  circular"  curve  C'  in  the  2 ' plane.  The  transforma- 
tion is  analogous  to  a Joukowski  transformation,  and  is  given  by 


(3-1) 


where 


Z = *c  * * /c 


(3-2a) 


(3-2b) 
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The  independent  variables  in  the  ^ ' plane  are  taken  as  ( ^ Oo  ). 
The  constant  is  defined  by 


*6  =.  c - 


fc  Ur 

" z 


~^Tu- 


CL 


where  and  ^Tc-  are  the  leading  and  trailing  edge  radii  of 

curvature  of  the  airfoil. 

Substituting  (3-2)  into  (3-1)  and  equating  real  and  imaginary 
parts  gives 


X • — Z.  b ol  ( o c^o-c  ^ 

, I 1 

/ - to  ^ 


(3-3a) 

(3-3b) 


This  transformation  is  unique  if  CO  is  restricted  to  Ot : 

Eq.  (3-3)  is  then  inverted  to  give  UJ  (■)<<,  } yc  ) and  ^/(Xc  yj,  ) • 
Combining  (3-3a)  and  (3-3b),  and  using  trigonometric  and  hyperbolic 
function  identities  yields 


r 


CO 


■t[» ;*x*(4r  J 


1 


(3-4a) 


and 


(3-4b) 


”here  = /-  />4 ) 


Xc 
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The  correct  sign  for  , and  the  correct  quadrant  for  LO  follow 
by  inspection  of  Eq.  (3-3). 

With  yy(yc  ^and  c/(*c,%)known,  the  first  transformation  is 
complete.  That  is,  for  each  airfoil  input  point  ( Xc , ) , a point 

( y ) (O  ) is  found.  The  resulting  points  ( 'f'  ( CO  ) are  considered 
to  form  a function  y>  =.  The  curve  C/  in  the  Z.1  plane  is 

then  given  by  % ' — b The  function  ^(Co)  is  gener- 

ally a slowly  varying  function  (this  is  equivalent  to  saying  is 
"nearly  circular")  and  is  periodic  with  period  Z-Y 

The  transformation  from  the  Z.'  to  the  ‘j  plane  is  now  con- 
sidered. Polar  coordinates  as  shown  in  Fig.  3 are  used,  where  (\(<p) 
are  taken  as  the  Independent  variables.  On  the  circle,  which  is  cen- 
tered at  the  origin,  X * XK  = constant,  and  the  radius  of  the 
circle  is  therefore  'L  « o c- 

The  transformation  used  is 

- ~ Y — - (3-5) 

/vi-  / 


as 


where  the  constants  ( ) are  real.  Since  zf/— 'v 

^ oO  , the  planes  coincide  at  oo  and  the  free  stream 

conditions  are  the  same  in  both.  Substituting  the  respective  polar 
coordinates  for  z'  and  ^ into  (3-5)  and  equating  real  and  im- 
aginary parts  of  the  resulting  expression  gives 


v-  y = 


A1  --  , 


and 


yvi  - i 


pr  v „ *o<Lsnip  & (p } 

' Y X ^ 

~gr~  ( ^ C*y4  /n.if  - O^..  x>  '.v\  yw  cp) 


~ (p 
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The  last  two  equations  represent  the  mapping  of  any  point  in  the  2. 
plane  to  the  corresponding  point  in  the  1 plane.  The  coefficients 
^ ^ ) are  determined  from  the  requirement  that  the  curve  vl/ 
maps  into  the  circle.  Therefore,  from  this  point  on  the  variables 
will  be  used  in  the  restricted  sense  that  ( oj)  define  the  curve 
C#  and  ( X , (p  ) define  the  circle.  Then,  since  X — X^on  the 
circle. 


no 

X*  + 21  ( A,*  C/S-4-  'll  Q + C //■ i/?  J (3-6a) 


and 


60  - cP  = ZT  ( CL  .CXrO  /jrucp  - -dtI  Jwn/ricp)  (3-6b) 


where 


— ^ X K 

lx*- 

3 Xfc 


(/■  p. 


Eq.  (3-6a)  is  the  standard  form  of  a Fourier  series,  and  it  follows 
that  the  coefficients  are  given  by 

XK  = fio  - J VtCP)^ 

o 

^.7T 

= ~4~  / XL£<L/)up  cj<p 

o 

2-(f 

~&ja  - y”  j C ) AUr  c/f 

O 
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Eqs.  (3-6a)  and  (3-6b)  give  and  CO  as  functions  of  (y  , 
and  the  resulting  points  { (p)  } represent  points  on  the 

curve  <2  • From  the  2 "2  * transformation  it  is  known  that  W (uj) 

also  represents  the  curve  Q1  ; therefore,  y and  (y  represent 
the  same  dependent  variable,  but  considered  as  a function  of  differ- 
ent independent  variables.  The  solution  for  the  transformation 

then  requires  a solution  for  the  functions  satisfying  (3-6a)  and 
(3-6b)  where  'A  (lx>)  is  a known  function. 

It  is  convenient  to  introduce  a variable  £ defined  by 

€ = CO  - cfi  (3-7) 


Like  (y  , £ may  be  considered  a function  of  either  Co  or  (f  » 

and  a similar  notation  is  employed,  i.e.,  £ — 6( <p)  = £ (tu) • 

It  follows  from  (3-6b)  that  £ is  periodic  with  period  2-71*  • 
Applying  (3-6b)  to  a specific  point  on  the  circle  ( p = cpf and  using 

the  definition  of  the  Fourier  coefficients  gives 

2-7/ 

“■(<£>')  - 2T  ( q)  vcf'1  J(p 

tT  1 Ao^rt<p' 

o 


Interchanging  the  order  of  summation  and  integration,  and  using  a 

trigonometric  identity  then  yields 

orr 


JOG 

^(<p)  z: 


vv  jA  (cp-(P‘)  c/(p 
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Introducing  the  rather  obscure  identity 


>-  .v(y-<p •*+[<**■')  z£  l 

V X i-  ^ ■ tp-cD‘ 


2- 

2- 


then  yields 


€^'}  = ^ d<p 

O 


/ r (U'^f (2.AJ+ ^ - 7 9 

" J ^ ;—  <p-P/  — JP  / 

o ^/Lasa.  ■ — j—  J 

The  first  integral  is  Independent  of  \s)  , while  the  second  is  zero 

in  the  limit.  The  equation  then  becomes 


Xff 

€(<p‘)=  r 00)  set  J<P  O-8) 

O 

The  integral  in  (3-8)  is  singular,  and  the  Cauchy  principle  value 
must  be  used.  The  integral  is  usually  called  Poisson's  integral. 

If  the  function  ^((p)  were  known,  (3-8)  would  be  a definite 
integral  for  £((p')-  However,  what  is  known  is  p (u))  , where 
p — CO  — £ , and  (3-8)  is  therefore  an  integral  equation  in  € . 

The  integral  equation  is  solved  by  the  method  of  successive  approxi- 
mations. Convergence  of  this  method  ultimately  depends  only  on  the 
function  p(.uS)a.nd  the  initial  approximation  for  the  function  £ . 
For  most  airfoil  shapes,  p (to)  is  a slowly  varying  function,  £ = o 
is  a good  initial  guess,  and  only  one  step  in  the  iteration  is  suffi- 
cient for  very  good  results.  In  this  case,  (3-8)  can  be  reduced  to 
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2-ir 

| p | - 

<£  6x>‘)  ~ 6,  c 00')  =•  ^ | c^t  • ^Clo)oJoj  (3-9) 

o 

F.q.  (3-9)  gives  £ (tu')as  a definite  integral  since  ^ (U>)  is  known. 
In  general  the  integral  must  be  evaluated  numerically  since  ^ (U>)is 
only  knovm  numerically.  This  is  done  by  prescribing  a value  of  C o' , 
performing  the  integration,  and  repeating  this  process  until  a suff-*- 
icient  number  of  points  are  obtained  to  define  the  £ (Oj)function. 
With  £ ( co)  known,  (f(io)  — do  — £(Ul)is  known,  ^ ( (p)  is 
therefore  known,  and  the  second  transformation  is  complete. 

The  velocity  \/  on  the  surface  of  the  airfoil  is  obtained  from 
thecomplex  potential  function.  Since 


ofUT(i) 

\yCi~  J 


airfoil 


it  follows  that 


^ ) 


airfoil 


V - 


(3-10) 


The  complex  potential  for  flow  over  a circle  of  radius  (X  is 


' ’P 


Z ,-,*<*• 


o^e.  \ 
J"  > 


■h 


lr_ 


ye 


-Zoi 


CL 


Therefore 

JtAr(r) 

dy 


The  lift  coefficient  is  related  to  the  circulation  p by  the 

Kutta-Joukowski  theorem,  which  gives 
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r = 


C \4, 


Then 


cJ'-U(J) 

~w 


circle 


— Z l/pt,  j +■  -7^— - 


Taking  differentials  of  H’and  J (with  ^ - \y.)  gives 

£1'  = (/+  f'EWu) 

2 ' \ / 


and 


o/  j"  • / o/  6 


r 


- i ('  - 


Then,  since  J rr  ^ (£  <3  ^ 


on  the  circle 


all.'  I c' 


be*h-  it! 


{z\'F+(Sr 


Differentiating  (3-1)  gives 


oi  i 


k2~ 


* - J-  ( 2'-  — ) 

7 - 2/  ^ 21  £*  J 


, uj  l CO 

Using  £=■  h 2 2 in  the  terms  in  parenthesis  then 


f 


_ l zM, 


I 


o/ir  ( 

7T*  I alr£o11 


(3-11) 


(3-12) 


yields 

(3-13) 


Eqs.  (3-10)  through  (3-13)  then  give  the  velocity  on  the  airfoil  sur 
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face  as 


14 


<oXk  1/2.-OA  (£0  ~C>0  +- 


C<C 


/ - 


ole  1 

&»  I 


(3-14) 


rVO 


'ytUwuk  V J 


Eq.  (3-14)  applies  to  airfoils  with  either  sharp  or  blunt  trailing 
edges.  For  an  airfoil  with  a sharp  trailing  edge,  ( O UO- O ) 
at  the  trailing  edge,  and  the  Kutta  condition  becomes 


( <p  - oO 


C^_ 

rr 


O 0l±  to  - o 


Since  (p  (^o)  =.  — £ (o)  > the  Kutta  condition  becomes 


<U  = (ec°)  •*•*0 


(3-15) 


<?> -it*  >•  e 


X, 


, &V"- 


( fe  to)  <r-  »--) 


The  pressure  coefficient  on  the  airfoil  surface  follows  from  the 
Bernoulli  equation,  and  is  given  by 


CT 


/ ~ 


/J4 

Wo 


o° 


) 


(3-16) 


The  potential  flow  is  then  known  since  [/y  and  on  the 

airfoil  surface  are  known.  A different  form  of  the  potential  flow 
is  required  as  input  to  the  boundary  layer  program,  however,  and  this 
is  considered  next.  Curvilinear  coordinates  shown  in  Fig.  1 are  used, 
and  ]/y  becomes  (XQ  ♦ The  curvilinear  coordinate  is  measured 


along  the  body  surface  from  the  front  stagnation  point,  and  it  is 
therefore  necessary  to  determine  the  arc  length  along  the  airfoil. 
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The  arc  length  /\^  between  two  given  points  ( )Cc 

( Xc  . Vr  . ) is 

) /CV+' 


« yc  • ) 

V > / 


and 


= 


rpm  jr‘  - r/w  * 

V */ 


The  first  form  of  this  equation  is  used  for  all  increments  except 

those  adjacent  to  the  leading  edge  (and  trailing  edge,  if  the  trail- 
ed Yc 


ing  edge  is  blunt),  where 


the  second  form  is  used,  but  not  directly  since 


is  usua- 


For  these  increments, 

o{  Vc 

7? 7c 

lly  known  as  a function  of  instead  of  . The  square  root 
term  is  expanded  to  give 

I /JXcs2-'  1 /^Xcx2-  / JXz 


Then 

yc 


Xc  • 

/ 


- 1 < *,•>  * t j'-tS 


V7»xJ 


Each  integral  is  evaluated  numerically.  With  the  arc  lengths  known, 
the  value  of  X is  calculated  for  each  output  point  (l.e.,  each 
point  at  which  the  velocity  \J . — CiQ  is  known). 


The  functions  V 
quired . 


and 


defined  by  (3-17)  are  also  re- 


*/c 


£*(c)-  f 


(3-17a) 
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✓*  /'J*  ' _ 
v c.  y 


ol 

^T) 


(3-17b) 


where  U*  / \J^> 

Both  "■  * and  are  evaluated  at  each  output  point  by  evalua- 

ting the  required  integrals  and  derivatives  numerically.  The  de- 
rivatives are  evaluated  using  least-square  parabolas  for  five  con- 
secutive nonequally  spaced  points. 
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4.  BOUNDARY  LAYER  ANALYSIS 


The  boundary  layer  properties  for  both  laminar  and  turbulent 
flow  are  calculated  using  the  Cebeci,  Smith  finite  difference  method. 
The  turbulent  Reynolds'  shear  stresses  are  eliminated  by  using  an 
eddy  viscosity  expression.  For  laminar  flow,  the  eddy  viscosity  is 
set  equal  to  zero,  and  the  equations  reduce  to  the  classical  laminar 
boundary  layer  equations.  Both  the  streamwise  and  normal  derivatives 
are  replaced  by  finite  difference  expressions,  and  the  resulting 
algebraic  equations  are  then  solved  by  a matrix  factorization  method. 
The  Cebeci,  Smith  method  has  been  extensively  developed  and  tested 
against  experiments  and  other  theories  and  has  been  found  to  be  acc- 
urate (Ref.  8).  The  method  also  predicts  either  laminar  or  turbulent 
separation.  The  Cebeci,  Smith  method  is  reviewed  herein  in  sufficient 
detail  so  that  the  computer  program  in  the  Appendix  can  be  understood. 

The  method  has  been  developed  for  either  two  dimensional  or  ax- 
isymmetric  flows  (with  or  without  transverse  curvature  effects). 

For  incompressible  turbulent  boundary  layers,  the  conservation  equa- 
tions are: 

Continuity: 


Momentum: 


[*6u 


2# 


- o 
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where  £ r>  for  two  dimensional  flows,  and  £ — / for  axis- 
ymmetric  flows.  The  quantity  /i(X}y)  Is  the  radius  from  the  axis 
of  symmetry  to  the  point  being  considered  in  the  boundary  layer. 

In  terms  of  the  body  geometry 

A (X,y)  = Cx)  + LO  y AUML  <p(x) 

where  Auj  (X)  is  t^ie  radius,  (p(X ) is  the  angular  inclln 

ation  of  the  body,  -=•  / if  transverse  curvature  effects  are 
Included,  and  £<J  — O if  transverse  curvature  is  neglected. 

The  boundary  conditions  considered  are 


(X,o)  = O 

U~  (X ) o)  n tfcj.  (x)  (suction  or  injection  permitted) 
6L  (X}*>)  = 

Defining  an  eddy  viscosity  £^by 


- (°  u!irl 


Av  e + 


]UL 

9/ 


and  using  the  Bernoulli  equation  for 
Momentum : 


9P 

3X 


then  yields 


94.  , 

/°&  w"  + /“</■ 


- A dt 


cku 

dy 


A/>' 

d/L 


(/.evfj- 


“O 


The  continuity  equation  is  identically  satisfied  by  introducing 
the  stream  function  ^ , defined  by 
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The  independent  variables  are  transformed  from  (X^y)  to  ^ 

where 

I)  00  J'  (4-la) 


PT 


/ 

/ /^/ 


(4-lb) 


and  a dimensionless 


¥ 


stream  function  -f  (^*  ^ — is 

^6  ' V ZS 

duced.  Then  =.  — , and  the  momentum  equation  becomes 


intro- 


+ rl£  f 2tf—  — - 2£  Z£1  - 


(4-2) 


where 


w - # ^ 


and 


, , ^ / A-  )‘tc 

(/+t)  = " /'*€" 

The  quantity  if'  ( ^ ' j is  the  transverse  curvature  term. 
The  boundary  conditions  transform  to 


/L  \2-£ 


pTf  /tO<L  (P  7 

, 


(4-3a) 
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fAilTv 

7W 


(4-3b) 


(4-3c) 


Eq.  (4-3a)  is  used  in  obtaining  (4-3b).  Eq.  (4-3b)  can  be  further 
simplified  by  solving  it  as  an  ordinary  differential  equation  in 


-f  o)  • Then  (4-3d)  replaces  (4-3b). 


= 


(4-3d) 


I ^Ah 


The  Cebecl,  Smith  eddy  viscosity  model  consists  of  an  inner  re- 
presentation near  the  wall  and  an  outer  representation  away 
from  the  wall.  The  inner  expression  is 


Z t- 


C 


*,  y 

■p' 


where  |(  is  a constant  and,  for  an  impermeable  wall*. 


Aa  = 


L6  v 


ft* 

r 


+ d ?e  Y 


dx.  r 


*See  footnote  on  next  page 
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The  outer  expression  is 


e?  = 


Kj.de  i u*)*1* 


where  /C  is  a constant.  The  boundary  layer  thickness  is  taken 


as 


a 


the  value  of  y where'-^  ■ 0.995.  With  this  definition  of  ^ , 

the  values  for  the  constants  are  taken*  as  ■ 0.4  and  K^m  0.0168. 
The  inner  and  outer  expressions  are  used  as  follows  (for  a given 
value  of  pi  ) : The  eddy  viscosity  is  taken  as  £ y for 

0 ~ * and  Is  taken  as  £■*  for  y 4 y £ § • The 

value  of  yc  is  the  value  of  y at  which  =-  £*.  The  eddy  vis- 
cosity £+*is  thus  continuous , but  in  general  its  derivative  has  a 


finite  discontinuity  at  y-=.  . In  terras  of  (^  ^)  * the  eddy  vis- 


cosity expressions  become 


<?+=  ^ci+ty 


/Tf 

^iytar 


X.  , n, 

Y Ifl  ♦ 


I - 


-EL  x j^L(fii./iY) 

« I {V?  '+4r  r 1 ' 


(4-4a) 


where 


7 


Y<$?>  = J 


7 

0*/ 


r//__  "5Z£>  n11  _ r11  / w \ 

1 -f  = yp.  1 -fur  ~ T (f,°) 


*More  recent  expressions  for  K^,  , and  A , are  given  in  Section  5. 
They  should  also  be  used  here,  but  became  available  too  late  to  make 
the  necessary  changes. 


A 


& 


l 
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<=t  - 


>u.K 


S! 


nr 


I + f.r 


£ (l+t)  d*\ 


(4-4b) 


For  two  dimensional  flow  or  axisymmetric  flow  without  transverse 
curvature,  ’t-O*  and  these  simplify  to 


<=  k'M  ^ irh 


I FT 


(4-5.) 


4 _ Kzji?  Vs  Z ~ fy) 

" x-AJi,  ,*r. 


(4-5b) 


In  the  present  application,  only  two  dimensional  airfoils  are  con- 
sidered, and  the  simpler  forms  (4-5)  are  used  (with  ^ = O also). 
It  is  convenient  to  nondimenslonallze  ^(X)  • This  is  done 

by  defining  ("$“)  as 

' X/c 

r(4)«  (4-« 


where  <£4,  = 


= <*//. 


Then 


* = r 


where  f [/^C  / ll  • The  momentum  equation  (4-2)  then  be- 


comes 
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Jy  [("if*  (HC1)  tI'O  * ['-  (l»rj 


2 


r H t 1 r*  I 2£  BI£  _ 2£  Zlf.  1 = 
* I If  •>  I 0|*  D*)  2^7  J 


where 


2-|*  olUl 
U*  M ?* 


and 


2<£  / <2  /I?7  AJ^lO 

<«,  - „ «»(*J  * 


The  eddy  viscosity  expressions  (4-5)  become 


€?=  Kf^jrSP^- 


r 


i - j-tfi 


\aJ 


«/t  ^ *v 


x,'.'  ( 

xt, 


< = K<7rSf^  W'  ^77^ 


" ^-(lf 


The  boundary  conditions  become 


O 


#• 


(4-7) 

(4-8) 


(4-9a) 

(4-9b) 

(4-10a) 
(4-1 Ob) 


(4-10c) 
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For  two  dimensional  flow  ( £ — o ) , 6^0  follows  from  (3-14) , 

and  ^and  £ reduce  to  the  definitions  given  in  (3-17). 

The  output  from  the  potential  flow  program  therefore  provides  all 
the  potential  flow  input  needed  for  the  boundary  layer  program. 

Eq.  (4-7)  subject  to  the  boundary  conditions  (4-10)  is  solved 
by  a finite  difference  procedure  (Ref.  8).  First,  the  ^derivat- 
ives are  replaced  by  finite  difference  expressions.  For  a general 
function  $ > a three  point  Lagrange  polynomial  gives 

££ 

d** 


= A,  }(&)  + «-n> 


J/V1 


where 


Jr*  - Vf* 

J>  /VX 


M-l 


** 

Os* 


= 


- u 


* _ <t?>( 

/V\  ^ yvi-  X 


) 


The  ^ * increments  do  not  have  to  be  equally  spaced.  For  a two  point 
formula,  (4-11)  still  holds,  with 
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I 


- 


Cl 


-/ 

V+-  _ 

'A! 


ft  3 = 0 

With  the  J-  derivatives  replaced  by  finite  difference  expressions, 
the  momentum  equation  (4-7)  becomes 


r ^ 


( 1+  ■?“  ] ' 


AiC 


where  primes  denote  ^ derivatives.  If  solutions  are  known  at 
“ ^ ( and  ^ t this  is  an  ordinary  differential 

equation  for  ( 7 ) 

Cebeci,  Smith  noted  that  round-off  errors  could  be  reduced  by 
introducing  a "translated  stream  function"  (p  defined  by 

P = f-  7 

Since  + A2  + “ 0,  the  momentum  equation  in  terms  of  (p is 
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(tit1)*  Q"1  ~ ( X<p'  + <p,Z) 

+-  V * (.V  * V) 

+ if*  [ *>*  ( A,  9+ Ai  <pm„  y a,  &..J 

~ ( f *4  l)  (A  tp'  ■*  Ax  (f'm -1  * Aj  (fijn-j.  )J~°  (4-12) 

where  Cp  = Eq.  (4-12)  will  be  solved  by  an  iterative  proced- 

ure. It  is  first  linearized  by  assuming  certain  terms  that  make  it 
nonlinear  are  known  from  a previous  iteration  (denoted  by  <pp)-  No 
terms  are  neglected,  however,  and  as  the  iteration  converges  the  sol- 
ution of  the  full  nonlinear  Eq.  (4-12)  is  obtained. 

The  linear  equation  is 

/i,  (l<p'  + + <p“  OPo*7!) 

5/*  £ ( A i (p  y A 3 @*1-2. ) 

/ (4-13) 

- OpJ+OiA'V  + **$.-.  1 Aj  &-*.)]  = 0 

Cebeci,  Smith  found  that  round-off  errors  could  be  further  reduced  by 
writing  (4-13)  in  terms  of  A (jp  — Cp  - The  equation  then  be- 
comes 

e,  A(p"'  + Apn  + Ez  A(f‘  t Eq  A <p  = €> 


(4-14) 
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where 

C,  = Oiil  ( It  6% 

^ * a 

e3  = -/>~  ( ?i  1 1)  - z$*  ( cpj+i)  fl, 
f,  = fi"  a, 

£>  = - ( f,  e,. #"  * e,  <Pi  * £V  <oJ 

+ [(  ft'  + 0 ( * * &'-.  +•  ^5  ) 

— &'  (Ax<Pa ,., 


With  values  of  known  from  the  previous  iteration,  Eq.  (4-14)  is  a 
linear  equation  for  £ <p  . It  is  solved  by  replacing  "/  derivatives 
by  finite  difference  expressions.  The  finite  difference  expressions 
are  obtained  from  the  five  point  Lagrange  polynomial 


3/.,  i L‘Cy) 


•/  L ' ('*)  ' y L ‘ (+)  9 • 


e+x 
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where 


LSV 


i?L  (i-™ 

<■¥-  f 


t+l 


K¥^  / 


Then 


f=  L'u  ?*•-«. f & ?«  +L\-v  + 4j,  ?«;. 

with  similar  results  for  the  higher  derivatives. 

A constant  step  size  in  the  ^ direction  can  generally  be 
used  for  laminar  flow.  For  turbulent  flow,  however,  a constant  step 
size  results  in  an  excessive  number  of  points.  This  is  because  tur- 
bulent boundary  layers  have  large  gradients  near  the  wall  (which  re- 
quire small  step  sizes) , but  may  be  many  times  thicker  than  laminar 
boundary  layers.  A simple  variable-step  grid  was  therefore  developed 
by  Cebeci,  Smith.  In  this  grid  the  ratio  of  lengths  of  any  two  ad- 
jacent intervals  is  constant.  Then 


$t‘  — 


K 


constant 


The  step  size  then  grows  progressively  larger  away  from  the  wall. 
Taking  - o and  ^ = P ^ , the  distance  to  a grid  point  is 
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7-,  = S', 

!M|  / 


AT1-/ 
K-  / 


The  variable  grid  involves  the  three  parameters  Hj  j K » and  the 
number  of  steps  . Table  1 in  Ref.  8 gives  suitable  combinations 

z 

of  these  parameters. 

Applying  the  finite  difference  formulas  to  this  grid  and  to 
(A-1A)  gives 


for  c = 3,  . . . 

N—  Z. 

where 

II 

r>  . -*■ 

M (a,  + K -KzdL,) 

+ M,( 

-§7  <a,  - k’")  ^-/c3af 

* = X\ 

' 6 (&-/  + K — Kz  - Kzat) 

+ L — 

e, 


K (k3ci,  -k\  ~Kz-  Ka*  - KQl,  * a,) 
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M 

r\ 

II 

7t)  S) 

^6(3,-/^''  ^(l,) 

* - f ,- 

^K'x.Ot'-CL,-,)  * £*,•,  K^,1! 

6K(i-K-Kd,)ii^  K1  (^a,  -a,  '0 

+ V,  ■ 

#,u  «'*} 

N-  = 


r T>  * 

62  / /1  *-;l 


= Br  ; 6 (a,  + k - )cij  + i£t  k (a,  - a,k-Kz) 


1. 


C3  ,r3 

J T, 


- /7/.,  ~ d, 


a<  = } ^t=  / f K •*  Kl  ; d3  = l+Z+lC'+k3 


"6,  = 


1 


1 A 


3*  = 


- 1 


x a,axK3tf^ 


3^'  ' 83  = 


« _ “J 

~ o j i/6 


d, 


1, a^K6  fa 
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A similar  procedure  for  the  point  /V-/  gives 


where  7?*-,  = "Bj-  [ (,  (<Xx+KX,  +KZ)  + 2.  ^ K‘ 

• v + a**  +a.,  + -|t-  4^v  K5  a,  -2^ y 


F = ^ 

t-N-,  £ 


-=T^  | (-U  + Kd,  ^ /EC 


•-/<-3/>  + i -c  C,  * • 


(4,at  f Aj/<-  u^K^  + a,  k%-  a.,^-  k9) 


t5  /'•> 


e/ 


z'J(':t/  4,)  + 


£/  By  /?. 


A'-y/ 


c - 3i 

Vi  a' 


/V-i  (''tl  ^ k3)+i  £i 


Cy  * * 


( 3cfu^ K & ) G\ 


aM  - §-  C- 


I = -!*-  f t + z </2 (a^-^K- I<V 

-Hj-i  ,cAA-/  C c/  " v 


^ ***■■*  2-  i 

/a'-v  J 


Mu-i  p 


^ ^/<  (d^K-^)  +1  ■^•Cy  * ' 
. (0.,-K0.,-K~)  ~ -ff  Cv 


^M-l  ~~ 


CjT 

^"/  ^A/-| 
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The  boundary  conditions  (4-10)  in  terms  of  Cp  become 


^u-  o/f* 

(^)4 

(4-17a) 

= -/ 

(4-17b) 

1U)®  ° 

(4-17c) 

Each  iteration  for  (p  satisfies  these  boundary  conditions,  so 
the  boundary  conditions  for  ^ (p  are 

AlPur  = O , 4<j C = O , &([>'s  = o 

In  the  finite  difference  scheme  the  subscript  ( ) is  identical 
to  the  subscript  ( ) hence  the  first  boundary  condition  provides 

Afy  - O (4-18a) 


Applying  the  finite  difference  procedure  to  the  second  boundary 
condition  and  using  Eq.  (4-18a)  gives 


&<9S  + f-  d<p,  &<PX  = ° 


(4-l8b) 


where 


G 


% 


/3y  ct/  <X_3_ 

Rz. 

Q-’X.  3 


AO 


Jl  - 

d,  ( x u j 

II 

"Br  $-/ 

Applying  the  finite  difference  procedure  to  the  last  boundary 
condition  gives 

4 ft-  + + Z.At-t  /M^<Pw.70(4'18c> 


where 


t~KJ  ~ 


= 

M*  = 

Kn 

^A,  = 

Brf 

KZ&,i 23 


- By  (k3&,cli_  * Kza,a.s  tKiias  + a,axa, ) 

With  (pt  — C>  known,  there  are  N-l  unknowns  A<pc'  } C — Z,  ■ . .{  N. 
Eq.  (A-15)  yields  N-A  equations  (one  for  each  point  C — 5^  • * 

Eqs.  (A-16),  (A-18b),  (A-18c)  each  yield  one  equation,  so  there  are 
N-l  equations. 

This  system  of  equations  can  be  rewritten  in  matrix  form 


[a]MJ  = M 


(A— 19) 
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where  [A]  Is  an  (N-l)  x (N-l)  matrix  consisting  of  the  /y  , 
Gc‘  } X'  f fA‘c  values,  [ N Is  an  (N-l)  x 1 matrix  consisting 
of  the  values,  and  is  an  (N-l)  x 1 matrix  consisting  of 

the  unknown  A Ql  values.  Eq.  (4-19)  Is  solved  by  a matrix  factor- 
ization technique.  Details  are  given  in  Ref.  8.  Once  (4-19)  is  sol- 
ved for  A , new  values  of  (p  are  calculated,  and  the  pro- 

cedure is  repeated  until  convergence  is  obtained.  The  convergence 
criteria  is  based  on  for  laminar  flow  and  on  (p ^ and 

for  turbulent  flow. 

Once  convergence  is  obtained  the  boundary  layer  properties  at 
the  station  v'+are  known.  The  velocity  profile  is  known  from 

0 /Y) 

—■  =.  -p1  — cp 1 +■  I . Other  parameters  (assuming  £ — O ) follow 

from 


cf  = 


KT 


fUe 


1 


/IF 


(4-20) 


s*  = - 
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(4-25) 


In  general  the  Integrals  are  evaluated  numerically.  The  value  of 

/ 

*94  is  the  value  of  for  which  j-  m .995  ( * -.005). 


Boundary  Layer  Solution 

The  solution  begins  at  the  stagnation  point  where  the  flow  is 


laminar.  A suitable 


^ grid 


for  laminar  flow  is  'ty  ■ 6 with 

/<*> 


^ ■ .025.  A solution  is  then  obtained  for  the  stagnation  point. 
Since  O at  the  stagnation  point,  (4-14)  shows  that  the  ^de- 

pendent terms  drop  out.  The  solution  is  not  completely  self -star ting 
however,  since  the  ^ terms  are  required.  A good  initial  guess  for 
y*0  was  found  to  be 

<p0'=  ‘i^'1  (°^h3) 


and 


<p.'  - 0 


(3&*/±6) 


// 


Values  of  (fi0  (p0"  ; (Po*  are  found  by  differentiation  and  in- 
tegration of  (p1 . The  iteration  continues  until  convergence  is  ob- 
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tained  (the  stagnation  point  solution  converges  to  the  Falkner-Skan 


solution).  With  the  stagnation  point  solution  known,  the  first  down- 
stream  point  ( J = ? ) is  calculated.  Two  point  £ formulas  are 

used  for  this  point,  while  three  point  formulas  are  used  for  all  sub- 
sequent points.  For  the  first  solution  at  a new  downstream  station, 
the  solution  for  the  previous  station  is  used  for  (pQ  ( (pj  (pj '! 

When  a solution  is  obtained  for  a given  point,  the  program  continues 
to  the  next  downstream  point.  This  is  continued  until  separation  or 
transition  occurs.  Separation  is  determined  by  — * O 

(i.e.,  / —>  O ) • A variety  of  transition  criteria  could  be 

0 7 lUr 

used.  Those  considered  for  this  analysis  are 

0/  Tq 


for 


O 


— 3 2-0  for  ~~  > O 

and  a criteria  due  to  Gaster  (Ref.  9).  Caster  obtained  a method  for 
determining  whether  laminar  separation  results  in  complete  separation 
or  a short  "laminar  separation  bubble"  followed  by  a reattached  tur- 
bulent boundary  layer.  If  the  laminar  boundary  layer  is  approaching 
separation  and  Gaster ’s  criteria  for  "short  bubble"  separation  is 
satisfied,  transition  to  turbulent  flow  is  assumed. 

When  transition  occurs,  the  ^ grid  is  redefined  for  turbulent 
flow.  Values  of  the  laminar  solution  at  the  previous  station  are 
obtained  at  the  turbulent  grid  points  by  interpolation.  The  solution 
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for  the  first  turbulent  station  is  obtained  using  two  point 
formulas,  while  subsequent  stations  use  three  point  formulas.  The 
turbulent  solutions  are  obtained  in  the  same  manner  as  the  laminar 

■f 

ones  except  that  is  finite  and  the  convergence  criteria  also  in- 
cludes S*.  After  the  first  iteration  (only)  at  a new  downsteam 
location  for  turbulent  flow,  the  old  and  new  solutions  are  averaged. 
Cebeci,  Smith  found  this  necessary  to  obtain  convergence.  The  tur- 
bulent calculations  continue  downstream  until  separation  occurs. 
Turbulent  separation  is  also  determined  by  the  condition  O . 

In  actual  calculations,  laminar  or  turbulent  separation  is  de- 

,11 

mined  when  (J^becomes  zero  or  negative  during  an  Iteration.  This 
Is  both  more  reliable  and  more  convenient  than  the  procedure  init- 
ially used  (Ref.  3).  This  procedure  cannot  be  used  with  the  Caster 
criteria,  however,  since  additional  Iterations  or  downstream  calcul- 
ations cannot  be  performed.  Therefore,  the  Gaster  criteria  is  used 

if  the  flow  is  laminar  and  ^ . 2.  . For  a boundary  layer  that 

II 

is  not  near  separation,  typical  values  of  0 are  > / . 

TUr 
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5.  WALL  JET  ANALYSIS 

The  turbulent  wall  jet  properties  are  calculated  using  a finite 
difference  method  based  on  the  Keller,  Cebeci  method  (Ref.  14).  Two 
dimensional  incompressible  flow  is  assumed,  and  a typical  wall  jet 
development  is  sketched  in  Fig.  2.  Wall  jet  calculations  are  a rea- 
sonably straightforward  application  of  boundary  layer  theory,  but 
there  are  two  difficulties.  The  first  is  that,  as  with  any  turbulent 
flow,  semi-empirical  expressions  are  required  to  evaluate  the  Rey- 
nolds stresses.  Such  expressions  are  available,  but  the  amount  of 
data  available  is  quite  limited  in  comparison  to  the  case  of  conven- 
tional boundary  layers,  and  the  available  data  are  not  always  in  ag- 
reement with  each  other.  This  situation  leads  to  uncertainties  in 
the  shear  stress  model,  but  an  apparently  reasonable  model  can  be 
formulated.  The  second  is  that  on  circulation  controlled  airfoils, 
the  wall  jet  develops  in  a region  with  a small  surface  radius  of 
curvature.  Boundary  layer  theory  can  account  for  this,  but  the 
boundary  layer  approximations  for  streamline  radii  of  curvature  begin 
to  break  down  over  the  latter  portion  of  the  jet  as  the  jet  thickens 
rapidly.  If  this  effect  is  not  adequately  accounted  for,  the  wall 
pressure  distribution,  and  therefore  the  separation  prediction,  will 
not  be  accurate.  An  exact  accounting  for  this  effect  is  difficult 
since  a required  correction  term  contains  second  order  streamwise  de- 
rivatives, which  do  not  fit  into  boundary  layer  theory. 

During  development  of  the  wall  jet  analysis,  calculations  were 
performed  for  a flow  denoted  as  Kind's  Flow  II  (Ref.  15).  This  flow 
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was  selected  as  a test  case  because  measured  wall  jet  profiles  along 
with  a measured  boundary  layer  profile  upstream  of  the  bloving  slot 
were  available  for  It.  The  effective  two  dimensional  angle  of  attack 
was  also  known.  This  flow  has  a moderate  momentum  blowing  coeffic- 
ient and  a relatively  small  trailing  edge  radius  of  curvature,  and 
should  be  a reasonable  test  case.  It  was  originally  intended  to  per- 
form calculations  for  other  flows  also,  including  those  of  Ref.  11, 
but  time  did  not  permit  this. 

Before  developing  the  wall  jet  equations,  the  boundary  condi- 
tions are  presented  since  they  are  simpler  but  illustrate  similar 
curvature  effects.  The  outer  boundary  conditions  require  that  the 
jet  flow  match  the  existing  external  flow.  When  the  surface  radius 
of  curvature  ^ is  small,  the  known  potential  flow  surface  values  do 
not  apply  directly  at  the  outer  edge  of  the  jet.  The  potential  flow 
is  unaffected  by  the  jet  according  to  boundary  layer  theory,  so  the 
following  analysis  yields  the  outer  edge  conditions.  The  irrota- 
tionality  and  continuity  equations  for  incompressible  two  dimensional 
flow  are,  respectively 


Integrating  the  irrotationality  equation  from  the  wall  to  the  edge  of 


O 


the  jet  gives 


^ - «e  = J fy 

o 

where  U g i8  the  known  potential  flow  velocity  at  the  surface  and 
Up  is  the  desired  outer  edge  velocity.  In  the  usual  analysis, 
as  in  Ref.  16  for  example.  (f  terms  are  neglected  since  U~ ~ 0 
at  the  wall  and  the  region  considered  is  thin.  In  this  case,  the 
continuity  equation  is  also  neglected,  and  the  outer  edge  velocity 
is 

(5-D 


Eq.  (5-1)  is  the  standard  curvature  correction  and  is  accurate 
within  a limited  distance  of  the  wall.  For  large  velocity  gradients 
and  relatively  thick  jets.  If  begins  to  become  significant  and  Eq. 
(5-1)  begins  to  become  inaccurate.  A higher  order  term  can  be  calcu- 
lated. Using  the  irrotationality  equation  with  - q gives  an 

expression  for  U.  , which  is  used  in  the  continuity  equation  to  ob- 
tain an  expression  for  If  . This  expression  for  U~  is  then  used  in 
the  irrotationality  equation  to  obtain  a new  expression  for  (JL 
For  « constant,  the  results  of  this  analysis  give  the  outer  edge 
irrotational  values  as 


4 e“  j ETaO" 


4- 


and 


ottu  o lliL 


(5-2) 


(5-3) 
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According  to  conventional  order  of  magnitude  estimates,  both  and 
the  second  term  In  the  brackets  In  Eq.  (5-2)  should  be  Insignificant. 


This  is  true  for  most  of  the  jet  region,  but  a combination  of  large 

'/U c 

r<y- 

ble  to  (J*  . The  term 


and  increasing  ^ result  in  (J^  eventually  becoming  compara- 

e 


MyC' 


in  Eq.  (5-2)  is  difficult  to  evaluate 


accurately  since  is  only  known  numerically.  It  is  calculated  by 


fitting  a least  - squares  straight  line  to  three  points  where 
is  known  and  using  the  slope  for  . 


c/ Uq 

My 


The  pressure  jpg.  at  the  outer  edge  of  the  jet  is  determined 


from 


(5-4) 


where  Is  total  pressure  at  the  outer  edge  of  the  jet  and 

jJ,o  and  (f^  are  known  from  Eqs.  (5-2)  and  (5-3).  For  irrota- 

tional  flow,  _7>  is  of  course  equal  to  the  known  free  stream 

’ -*-e  stag  n 

total  pressure.  However,  the  remains  of  the  upstream  boundary  layer 
can  result  in  a total  pressure  deficit,  and  this  can  be  accounted  for 
by  allowing  ^ to  vary  from  one  streamline  to  another.  This 

variation  is  relatively  minor,  and  will  be  discussed  shortly. 

The  above  analysis  is  used  for  the  outer  edge  conditions  for 
wall  jets  without  a velocity  minimum,  i.e. , after  the  minimum  has  been 
essentially  entrained.  Theoretically,  it  could  also  be  used  at  the 
outer  edge  of  jets  with  a velocity  minimum,  and  this  would  be  some- 
what simpler  since  ~~P.  would  be  equal  to  the  free  stream  value. 

For  typical  circulation  controlled  airfoils,  however,  there  are  prac- 
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tical  difficulties  that  suggest  a different  approach.  The  difficul- 
ties arise  because  the  upstream  boundary  layer  is  normally  several 
times  thicker  than  the  blowing  slot,  which  means  that  the  curvature 
corrections  would  have  to  be  applied  much  further  from  the  wall. 

This  would  be  no  problem  if  the  surface  radius  of  curvature  did  not 
change  abruptly. 

However,  most  circulation  controlled  airfoils  have  a discontin- 
uous radius  of  curvature  at  the  blowing  slot.  Applying  Eq.  (5-1)  or 
(5-2)  at  such  a point  would  result  in  a sudden  decrease  in  (by 
about  30Z  for  Kind's  Flow  II,  for  example).  The  radius  of  curvature 
could  be  faired  in,  but  there  is  little  information  to  suggest  the 
proper  fairing,  so  the  result  would  be  somewhat  arbitrary.  For  this 
reason,  and  others  discussed  below,  the  following  procedure  is  used. 

The  initial  wall  jet  calculations  include  the  velocity  minimum, 
and  the  outer  edge  conditions  apply  at  a point  corresponding  to  the 
outer  edge  of  the  upstream  boundary  layer.  The  outer  edge  conditions 
are  taken  as  the  uncorrected  potential  flow  surface  values 

£4=  Ue  , Je=  ?e  (5-5) 

To  compensate  for  the  lack  of  a curvature  correction  for  the  pressure, 
the  normal  pressure  gradient  in  the  wall  jet  equations  is  set  equal  to 
zero  for  the  portion  of  the  flow  above  the  velocity  minimum.  The  nor- 
mal pressure  gradient  is  actually  set  equal  to  zero  slightly  above 
the  minimum  Instead  of  at  the  minimum,  but  the  difference  is  not 
great.  Then  for  the  initial  wall  jet  development,  the  pressure  at 
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the  velocity  minimum  is  essentially  the  potential  flow  surface  pre- 
ssure J>&  , and  the  velocity  CL^ , w is  determined  during  the  wall 

jet  solution.  The  minimum  velocity  is  zero  at  the  slot,  and 

the  ratio  increases  as  the  jet  procedes  downstream. 

This  ratio  increases  rapidly  at  first,  and  more  slowly  downstream. 

As  the  jet  proceeds  downstream,  the  point  of  minimum  velocity  moves 
further  away  from  the  wall.  At  each  downstream  location,  the  known 
value  of  ^/h/aj  Is  compared  to  a velocity  calculated  by  applying  Eq. 
(5-2)  at  the  location  of  the  minimum  velocity.  Initally  is 

less  than  this  velocity,  but  eventually  the  two  values  become  approx- 
imately equal.  When  the  two  velocities  become  approximately  equal, 
wall  jet  calculations  switch  over  to  the  case  of  a wall  jet  without 
a minimum  velocity,  i.e. , the  flow  above  the  minimum  is  dropped  from 
wall  jet  calculations.  At  the  switching  point,  the  known  velocity 
profile  for  the  portion  of  the  jet  above  the  minimum  is  used  to  tab- 
ulate streamline  constants  and  their  corresponding  total  pressures. 
These  tabulated  values  are  then  used  in  Eq.  (5-4)  to  determine  g 

for  the  downstream  calculations. 

This  procedure  allows  a smooth  transition  from  profiles  with  a 
minimum  velocity  to  those  without,  and  allows  the  remains  of  the  up- 
stream boundary  layer  to  be  accounted  for  approximately.  By  the  time 
the  switch  is  made  to  profiles  without  a minimum,  the  velocity  above 
the  minimum  is  nearly  constant  and  the  shear  stress  is  thus  essen- 
tially zero.  The  total  pressure  is  therefore  nearly  constant  along 
a streamline,  although  it  varies  somewhat  from  one  streamline  to  an- 
other. A more  straightforward  analysis  would  be  to  continue  to  per- 


to  the 


form  wall  jet  calculations  out  to  the  irrotational  flow,  i.e., 
streamline  corresponding  to  the  outer  edge  of  the  upstream  boundary 
layer.  This  would  require  that  calculations  be  performed  out  to  a 

C/t. 

point  further  from  the  wall,  and  would  not  be  difficult  if 


enough  for  the  boundary  layer  approximations  for  streamline  radii  of 
curvature  to  begin  to  break  down.  As  this  occurs,  uncertainties  in 
the  outer  edge  conditions  and  in  the  wall  jet  equations  Increase, 
and  these  uncertainties  increase  as  the  distance  from  the  wall  in- 
creases. This  effect  is  severe  enough  to  Justify  calculating  only 
as  far  from  the  wall  as  is  necessary.  Therefore,  the  present  proce- 
dure is  adopted. 

Fortunately,  the  inner  boundary  conditions  are  simpler  to  use 
and  to  describe.  At  the  wall,  the  no  slip  condition  and  an  imperm- 
eable surface  result  in 


The  boundary  layer  equations  for  the  time-mean  values  for  two 
dimensional  incompressible  turbulent  flow  with  a small  surface  radius 
of  curvature  are 


were  sufficiently  small.  However, 


eventually  becomes  large 


u - o , (/■  - o 


(5-6) 


Continuity 
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where  ^ = -l^y®  — — — \^7y  &y  — the  shear  stress. 
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This  form  of  the  momentum  y equation  results  from  the  customary 
boundary  layer  approximations,  and  is  equivalent  to  assuming  a stream- 
line radius  of  curvature  equal  to  ( R+  y ).  This  is  adequate  over 
most  of  the  jet  region,  but  begins  to  break  down  over  the  latter  por- 
tion where  the  jet  thickens  rapidly  due  to  the  strong  adverse  press- 
ure gradient.  Experimental  results  indicate  that  the  pressure  change 
across  the  jet  tends  to  zero  as  separation  is  approached.  However, 
the  above  equation  predicts  a finite  pressure  change  across  the  jet. 
This  means  that  if  the  pressure  at  the  outer  edge  of  the  jet  is  corr- 
ect and  the  above  equation  is  used,  the  pressure  gradient  at  the  wall 
will  be  too  mild,  which  means  that  the  predicted  separation  point 
will  be  too  optimistic,  i.e.,  too  far  downstream.  The  effect  on  cal- 
culated results  can  be  severe  enough  that  separation  is  not  predicted 
at  all,  even  though  upstream  calculated  values  are  in  good  agreement 
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with  experiment.  It  is  necessary  to  obtain  a more  suitable  form  of 
the  momentum  y equation. 

Observation  of  various  calculated  results  indicates  that  the 
pressure  change  across  the  jet  must  be  accurately  predicted,  but 
that  the  form  of  the  momentum  y equation  used  to  do  this  is  not 
particularly  important.  In  other  words,  the  pressure  distribution  at 
the  wall  is  much  more  important  than  the  details  of  the  streamline 
curvature.  However,  in  order  to  predict  the  wall  pressure  accurat- 
ely without  resorting  to  experimental  data  for  each  case,  it  is 
necessary  to  consider  the  streamline  curvature. 

The  momentum  y equation  with  all  convection  terms  retained 
but  all  stress  terms  neglected  is 


_L  af 
t ly 


u'2'  lr2!£  „ 2K  _ 

“ dy  + ■£+/  ^ 2* 


According  to  boundary  layer  order  of  magnitude  estimates,  the  last 
two  terms  are  two  orders  of  magnitude  less  than  the  first  two  terms. 
Also,  some  of  the  neglected  stress  terms  are  the  same  order  of  mag- 
nitude as  the  last  two  terms,  so  this  form  of  the  equation  is  not  en- 
tirely consistent.  However,  the  effect  that  needs  to  be  accounted 
for  in  calculations  seems  to  be  due  to  streamline  curvature  rather 
than  additional  stress  terms,  so  the  above  equation  should  be  an  ad- 

7j- 

equate  starting  point.  The  (j~  — — term  could  be  included  in  exact 

? / 

form  without  difficulty.  However,  including  it  but  not  the  - 

a- 

dir 

term  results  in  a negligible  improvement.  The  term  is  more 

Q A 

difficult  since  it  results  in  a second  order  streamwise  derivative. 

It  is  not  difficult  to  write  a finite  difference  expression  for  the 
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second  derivative,  although  the  expression  must  be  centered  one  grid 
point  upstream,  but  the  resulting  calculations  become  unstable  at  an 
incredible  rate.  An  exact  curvature  correction  is  therefore  unten- 
able. 


The  second  order  streamwlse  partial  derivative  can  be  reduced 
to  ordinary  second  derivatives  of  profile  parameters  if  an  approxi- 
mate velocity  profile  is  used,  or  if  local  similarity  is  assumed.  If 
an  approximate  velocity  profile  is  used,  the  result  would  be  equival- 
ent to  that  of  Kind  (Ref.  17),  although  the  procedure  is  somewhat 
different.  However,  in  using  this  procedure  in  conjunction  with  an 
integral  method.  Kind  (Ref.  18)  found  that  stability  problems  result- 
ed as  separation  was  approached,  and  the  method  had  to  be  abandoned 
over  the  latter  portion  of  the  flow.  A closely  related  but  somewhat 
simpler  method  is  to  assume  local  similarity  of  the  velocity  profiles. 
This  is  equivalent  to  using  an  approximate  velocity  profile,  but  ig- 
noring some  of  the  resulting  terms. 

U~ 

Including  the  r=r-r  term  in  exact  or  quasi-exact  form  tends  to 
J A 

produce  stability  problems,  while  ignoring  the  term  tends  to  produce 
an  insignificant  correction  and  therefore  inaccurate  wall  pressure 
calculations.  An  approximate  but  plausible  treatment  of  the  term 
thus  seems  in  order.  Examining  Kind's  measured  velocity  profiles, 
which  seem  reasonably  typical,  suggests  a possibility.  For  y h.  y**, 
— — has  a large  negative  value,  while  for  V ? V , may  be 

positive  or  negative,  but  tends  to  have  a much  smaller  magnitude.  It 


follows  from  the  continuity  equation  that  U~  increases  fairly  rap- 
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idly  from  the  wall  to  y^  , then  changes  more  slowly.  Therefore, 
should  be  a reasonably  representative  value  for  the  outer  portion  of 


the  jet,  where  most  of  the  pressure  change  across  the  jet  occurs. 


If  so, 


cW/fA 


rli- 


- should  provide  a reasonable  average  value  of 


Jjf 

2* 


, and 


a workable  form  of  the  momentum  y equation  should  be  obtainable. 
Using  the  definition 


along  with  Leibnitz's  rule  and  the  continuity  equation  gives 

/2.  + _ ot  f h ^ y 

— ^ ~ Tty.  » 


Si 

Since  the  definition  of  5 is  not  modified  for  curvature,  the  — ^ — 
term  should  probably  be  set  equal  to  one,  but  the  difference  tends 
to  be  minor.  From  this, 

Hx  ~ alx1  ' ~3TyT  A"  - -TfJ  ~otx 

Dropping  the  second  derivatives  gives 


ol  (f/fA  _ ofCljv*  oW ** 

ol%  ~ ~dHC  -^c 

Using  this  as  an  average  value  of  _2iL  in  the  momentum  y equation 
produced  good  results  for  a test  run,  but  the  procedure  is  not  en- 
tirely satisfactory  as  it  stands.  The  reason  is  that  calculations 
can  tend  to  ignore  the  term.  If  Cl  happens  to  be  too  large  during 
an  iteration,  is  too  small,  and  tends  to  be  too  small 

“dTT 
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also.  Thus  the  correction  term  tends  to  be  too  small.  If  CL/y*  is 
too  large,  the  pressure  change  across  the  jet  tends  to  be  too  large, 
and  the  correction  term  being  too  small  tends  to  reinforce  the  pro- 
blem. 

diir 

Including  second  derivatives  would  tend  to  give  a larger  , 

o(% 

but  can  cause  stability  problems  and  probably  should  be  used  only  as 
a last  resort.  At  least  part  of  the  stability  problem  encountered 
with  second  derivatives  was  due  to  a convergence  problem  (described 
later)  which  tended  to  cause  oscillations.  This  problem  was  later 
corrected,  but  by  then  approaching  deadlines  did  not  permit  re- 
evaluation  of  various  procedures. 

The  full  momentum  y equation,  with  the  aid  of  the  continuity 
equation,  can  be  written  as 


r dy 


ur_ 

£ty 


Dropping  the  if  ^ term  and  using  an  average  value  for  (~£L  / 


gives 


-L  2Z  _ _*L  Ax, 7 = o 


where  -f  (X)  - TiL  [_  jy,  ( u ')L,  As  discussed  above,  there  is  a 

"Vi  VQ 

reasonable  hope  of  finding  a suitable  expression  for  j (%)  based 
on  average  values.  Because  of  time  limitations,  a simple  order  of 
magnitude  value  determined  from  the  potential  flow  values  is  used. 


Using 
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LT  ~ X dUt 

~ ~ gives 

7)  (tr  \ ^ X , _/_  /V^c  *• 

^ U£\cJJ) 

An  average  value  Is  taken  as 


which  gives 


■f  (x) 


2A 

X 


( °ld<L  '1^’ 

\ dx  J 


This  approximate  procedure  at  least  allows  useful  calculations  to  be 
performed  without  the  use  of  experimental  data.  Also,  as  mentioned 
previously,  the  normal  pressure  gradient  is  set  equal  to  zero  above 
the  velocity  minimum.  This  can  be  accomplished  by  setting 
above  the  velocity  minimum,  and  doing  so  does  not  result  in  any  di- 
fficulty in  the  calculations. 

The  continuity  equation  is  Identically  satisfied  by  introducing 
stream  function  ^ defined  by 


a 
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For  the  wall  jet,  there  is  no  particular  advantage  in  Introducing  a 
transformation  similar  to  that  used  for  the  boundary  layer  calcula- 
tions - Eq.  (4-1).  The  independent  variables  are  therefore  retained 
as  ( X 1 y ) * The  dependent  and  independent  variables  are  nondimensi- 
onalized  and  made  of  0 (1)  by  use  of  the  following  definitions 


y 

7=  -r 


e - c-  _ V 

he  — / ) F — 


T pv£ 


/ -e  - 3L 

^ ~ ^-  / F ) <C  - slot  thickness  (5-7) 

Velocities  nondimens ionalized  by  V^,,  and  lengths  nondimension- 
alized  by  C , are  denoted  by  an  asterisk,  such  as  (Je  = Ug  j and 

£*=-<l 1/C 


With  these  definitions, 
and  the  equations  become 
Momentum  x 


//*  — n*  EE.  tr*~  - — 
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(5-8) 


Momentum  y 


o ? /-/-<?  7 e 
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The  inner  boundary  conditions  (5-6)  become 


r=  © 


ai  ^ - o 


(5-10) 


The  condition  if—  O at  y — O results  in  the  condition  U £ t~  m 

constant  at  ^ =.  O , and  the  first  condition  in  Eq.  (5-10)  there- 
fore results  by  setting  the  constant  equal  to  zero.  The  outer 
boundary  conditions  are  asymptotic  but  in  practice  are  applied  at  a 
finite  location  y — A ( ^ ~ become 


BF 

'd'tj 


Qpe  ai  'I’  ld 


(5-11) 


Eqa.  (5-2)  to  (5-4)  give 


R*  ot7t& 
~ dp 


(5-12a) 


C2 

ZstAG 


(5-12b) 


For  the  case  where  Eq.  (5-5)  is  used,  the  values  are 

2. 

U * — Lie*  ) Cj>^  ~ C.7»  = /-  U*  (5-13) 

If  an  eddy  viscosity  model  is  used,  the  shear  stress  is 


= -v  F (i+e4) 
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i 


or 

/ v / ~b(A.  U.  \ 

r=-  v?  (ne+)  ( -57  - 

where  £+i»  the  dimensionless  eddy  viscosity  coefficient.  As  de- 
veloped shortly,  an  eddy  viscosity  model  is  used,  but  with  a modifi- 
cation to  allow  the  shear  stress  to  be  negative  at  the  velocity  max- 
imum. The  shear  stress  is  therefore  taken  as 

where  ^ is  related  to  the  mixing  length.  In  dimensionless  form 
this  becomes 


l _ Y*(j>L  , j] 

/+n  37  J J 

17  (5-14) 


Different  eddy  viscosity  expressions  are  required  for  different 
portions  of  the  wall  jet  profile.  These  expressions  are  developed 
below.  For  some  of  the  development,  it  is  necessary,  or  at  least  con- 
venient, to  know  the  approximate  velocity  profile  shape,  so  this  is 
considered  first.  Sketches  of  profiles  showing  notation  are  given  in 
Figs.  4 and  5.  One  of  the  better  integral  methods  for  calculating 
wall  jet  development  is  that  of  Gartshore  and  Newman  (Ref.  19),  which 
is  also  used  with  modifications  for  curvature  by  Kind  (Ref.  18).  The 
velocity  profiles  used  in  these  analyses  are 

/ y \ ^ 

= O'- 


(5-15.) 
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and 


u.  - CL%+  IJ^L  <£  ' L*  J 


(/^Y^rs)  <5-l5b> 


where  U^—  LIq^ot  profiles  without  a minimum  velocity 

and  U.^  — <UM/JLot  profiles  with  a minimum  velocity.  The  constant 
h — . The  exponent  //  varies  from  about  1/11  to  1/2,  with 

N ••  1/2  typically  taken  as  a separation  criteria.  Eq.  (5-15b)  pre- 
dicts an  asymptotic  velocity  decay,  but  experiment  (Ref.  20)  and 
calculated  values  when  intermittency  is  included  indicate  that  O.— 
at  about  (/'Xw)  / ° ~ 2-j  • The  above  expressions  are  used 

anytime  approximate  velocity  profiles  are  required. 

The  eddy  viscosity  in  the  region  y^y^s  determined  from 
Dvorak’s  model  (Ref.  21).  The  data  used  for  this  region  are  from 
profiles  without  a velocity  minimum,  but  the  results  (excluding  the 
intermittency  profile)  are  assumed  to  apply  for  profiles  with  a min- 
imum also.  Dvorak's  model  is  based  on  measurements  by  Wygnanski  and 
Fiedler  for  a variety  of  turbulent  shear  flows.  An  eddy  Reynolds 
number  Is  defined  as  in  Eq.  (5-16)  and  found  to  have  a constant 

value  of  about  15. 


In  general,  JJ ^ is  a velocity  deficit  and  u~  is  the  standard  dev- 
iation of  the  intermittency  profile,  and  thus  is  a measure  of  the 
size  of  the  large  eddies  (Ref.  22).  For  a wall  jet,  fj^  = is 

known,  and  if  a suitable  expression  for  is  available  Eq.  (5-16) 


LU  <r 

is  £+ 
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(S 


(5-16) 
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provides  an  eddy  viscosity  model.  Dvorak  gives  results  for  based 
on  data  by  Gartshore,  but  as  discussed  below  the  results  do  not  seem 
to  be  general  since  they  essentially  prescribe  a wall  jet  thickness 
that  is  not  always  in  agreement  with  experiment.  The  reference 
Dvorak  cites  (Gartshore' s dissertation)  was  not  available  to  the 
author,  but  another  paper  by  Gartshore  (Ref.  22)  contains  apparently 
the  same  data.  In  Ref.  22,  Gartshore  uses  measured  intermlttency 
profiles  for  a wall  jet  in  still  air  and  three  wall  jet  flows  in 
specialized  pressure  gradients  for  which  self-similar  solutions  exist 
to  provide  an  experimental  verification  of  Townsend's  large  eddy  eq- 
uilibrium hypothesis  (Ref.  23).  The  data  in  Ref.  22  can  also  be 
used  to  check  Eq.  (5-16),  and  the  results  are  in  good  agreement,  giv- 
ing a value  of  P that  varies  between  14.5  and  17.  The  data  for  <J~~ 
eT 

can  be  put  into  different  forms,  and  the  method  used  here  is  to  de- 
termine U- />Loas  a function  of  U8  / . The  following  curves  are 

used  to  do  this 
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The  mean  location  of  the  intermlttency  profile  y is  also  determined 
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by  curve  fitting  the  data  of  Ref.  22  to  obtain 
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For  a given  velocity  profile,  and  /_o  are  known.  Eq. 

(5-17)  then  provides  Q~  and  Eq.  (5-16)  then  provides  . The 

value  of  y is  determined  from  Eq.  (5-18).  The  value  of  is 
multiplied  by  a curvature  correction  /^described  later.  For  a wall 
jet  profile  without  a velocity  minimum,  the  eddy  viscosity  for 
is  taken  as 


^ = *«ei  r 


(5-19) 


where 


(5-20) 


is  the  intermittency  function. 

For  a profile  with  a minimum  velocity,  the  eddy  viscosity  dis- 
tribution for  y is  taken  as  follows. 

^ + — 1/ 


(&*  * y ± ^ ) (5-2D 
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(5-22) 


where  A » —*12.7  A and  y is  given  by  Eq.  (5-20) 
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(5-23) 


The  value  in  Eq.  (5-21)  is  the  same  as  in  Eq.  (5-19)  except  that  the 
intermittency  function  is  excluded.  Values  for  Eq.  (5-21)  are  de- 
termined the  same  as  for  Eq.  (5-19).  The  region  in  which  each  of  the 
above  expressions  is  used  is  somewhat  arbitrary  but  seems  reasonable. 
Eq.  (5-15b)  shows  that  the  maximum  negative  shear  stress  occurs  at 
about  )^a+-.«5‘/LoS0  the  region  for  Eq.  (5-21)  extends  somewhat  above 
that  point.  This  allows  a smooth  switch  to  profiles  without  a min- 
imum. The  region  in  which  Eq.  (5-22)  is  used  corresponds  approxima- 
tely with  the  portion  of  the  flow  above  the  velocity  minimum.  In  Eq. 
(5-22),  ^ogi.  “nxl®0®  value  in  the  upstream  boundary  layer 

at  the  slot,  and  is  known  from  the  known  boundary  layer  properties. 
The  approximate  values  for  y and  in  Eq.  (5-22)  are  given  by 
Dvorak  as  asymptotic  values  for  a boundary  layer  with  a large  shape 
factor.  They  also  give  reasonable  agreement  with  the  intermittency 
distribution  used  in  the  upstream  boundary  layer  calculations.  Eq. 
(5-23)  is  an  arbitrary  fairing  to  connect  the  other  two  expressions. 
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For  wall  jets  on  convex  surfaces,  the  radial  gradient  of  momen- 
tum in  the  outer  region  is  negative,  and  this  tends  to  produce  addi- 
tional turbulence  in  the  outer  region.  The  curvature  correction  j£c 
is  used  to  account  for  this.  Kind  (Ref.  18)  derives  a curvature 
correction  based  on  data  for  the  growth  of  a curved  wall  jet  in 
still  air.  Kind  uses  the  Cartshore  and  Newman  (Ref.  19)  expression 
for  eddy  viscosity,  which  is  essentially  the  same  as  Eq.  (5-16)  ex- 
cept that  *s  determined  as  the  sum  of  two  terms.  Kind  corrects 

both  terms  for  curvature,  and  an  identical  procedure  does  not  apply 
with  Eq.  (5-16).  However,  the  major  portion  of  this  correction  is 
retained  by  taking 


Kc  = [it  1.19  - ,.3 V ( ] 

*-  7S.  ' T<  7 J 


[it  /•  3¥  '■  [l  t . 3 J 1 


(5-24) 


This  correction  can  be  quite  significant  since  it  can  easily  increase 
the  eddy  viscosity  by  a factor  of  two  or  more. 

Eqs.  (5-17)  and  (5-18)  are  used  in  place  of  those  given  by 
Dvorak  since  Dvorak's  expressions,  when  used  with  Eq.  (5-20),  result 
in  the  eddy  viscosity  becoming  zero  at  about  the  middle  of  Kind's 
measured  wall  jet  profiles.  A similar  difficulty  also  arises  for 
Jones'  measured  measured  profiles  (Ref.  24).  This  apparently  is 
related  to  the  fact  that  experiments  often  show  different  wall  jet 
growth  rates,  especially  near  the  blowing  slot.  Use  of  Eqs.  (5-17) 
and  (5-18)  allows  the  Jet  thickness  to  be  determined  during  the  cal- 
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culatlons,  and  calculated  results  indicate  that  the  thickness  of  the 
outer  portion  of  the  jet  is  essentially  determined  by  the  maximum 
value  of  the  eddy  viscosity.  As  would  be  expected,  Increasing  the 
eddy  viscosity  increases  the  thickness  of  the  jet.  Dvorak  uses  a 
diffusion  equation  to  allow  the  actual  eddy  viscosity  to  lag  behind 
, and  this  procedure  was  used  initially.  However,  for 
Kind's  Flow  II,  this  procedure  seemed  to  limit  the  eddy  viscosity 
to  too  small  a value,  while  use  of  resulted  in  good  agreement 

with  experiment.  Therefore,  the  local  value  Is  used  in  the 

calculations. 

The  inner  region  ( of  a wall  jet  profile  resembles  a 

conventional  turbulent  boundary  layer,  and  the  conventional  eddy  vis- 
cosity model  for  a boundary  layer  thus  seems  to  be  a reasonable  start- 
ing point  for  the  wall  jet  formulation.  In  general,  however,  such  a 
model  will  not  match  the  eddy  viscosity  in  the  outer  region  without 
the  aid  of  a somewhat  arbitrary  fairing.  The  conventional  boundary 
layer  model  near  the  wall  is 


(5-25) 


where 


Kz-b> 


while  away  from  the  wall  the  conventional  model  is 


o 
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(5-26) 
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The  matching  condition  for  the  two  expressions  is  that  the  eddy  vis- 
cosity be  continuous.  Common  values  for  the  constants  are  — .*/  * 
£ = .c/C  8 and  ^3=2^*  For  low  Reynolds  numbers,  Cebeci 

(Refs.  25  and  26)  recommends  the  values 
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(5-28) 

(5-29) 


Eqs.  (5-27)  and  (5-28)  apply  for  3oo  , while  Eq.  (5-29) 

applies  for  1^2^  H7.£~%  where 


o 


The  above  values  are  used  here,  and  if  Is  less  than  the 

given  limiting  values  the  expressions  are  evaluated  for  the  limiting 
values.  The  — term  in  the  expression  for  /Ao  Is  evaluated  in  terms 

f 

of  and  the  pressure  gradient  by  the  usual  procedure.  Near  the 

1° 

wall  the  convective  terms  are  neglected,  and  the  momentum  ^ eq- 


uation is 
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Evaluating  the  pressure  gradient  at  the  vail  and  integrating  gives 


This  expression  could  be  used  as  is,  but  for  conventional  boundary 
layers  Cebeci  (Ref.  28)  recommends  evaluating  it  at 
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The  expression  was  used  with  y as  a variable  and  also  with  y ev- 
aluated at  various  fixed  locations,  and  the  difference  in  the  solu- 
tion was  insignificant.  The  form  used  gives 


Aa  = 


*3^ 


(*♦%>]*  (5-30b) 


were  is  given  by  Eq.  (5-30a) . For  a convex  radius  of  curvature, 

both  the  inner  and  outer  eddy  viscosity  expressions  can  be  corrected 

2 

by  multiplying  by  $ , where  (Ref.  27) 
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Gartshore  and  Newman  (Ref.  19)  use  a correction  to  the  outer  eddy 
viscosity  to  allow  for  the  turbulence  in  the  outer  portion  of  the 
wall  jet.  This  correction  is  equivalent  to  multiplying  Eq.  (5-26) 
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by  the  factor  / , where 


Gartshore  and  Newman  use  this  expression  in  conjunction  with  an  in- 
tegral method  to  determine  the  shear  stress  at  y = , and  how 

it  applies  in  the  present  formulation  should  probably  depend  on 
where  the  fairing  to  join  the  outer  region  begins. 

The  eddy  viscosity  distribution  is  taken  as  follows.  Near  the 
wall,  Eqs.  (5-25,  27,  28,  30)  are  used  to  evaluate  . This  ex- 
pression is  used  until  , where  is  determined  from  Eqs. 

(5-26,  29).  Matching  occurs  at  about  y — .1)/^.  From  the  matching 
point  to  , the  eddy  viscosity  is  determined  from  the  somewhat 

arbitrary  fairing 
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A typical  eddy  viscosity  distribution  as  determined  by  the 
above  procedure  is  shown  in  Fig.  4,  and  is  drawn  to  scale  for  one  of 
the  calculated  velocity  profiles.  A lack  of  detailed  information  for 
the  inner  region  along  with  the  difference  in  magnitude  between  € ^ 
and  leads  to  some  uncertainty  in  the  formulation. 

Several  variations  for  and  the  fairing  were  used  during 

development  of  the  analysis.  These  variations  were  made  while  att- 
empting to  sort  out  development  difficulties  that  were  related  to  a 
convergence  anomaly  and,  eventually,  the  inaccurate  pressure  gradient 
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over  the  latter  portion  of  the  jet.  Most  of  the  variations  were 
made  before  these  difficulties  were  sorted  out,  so  conclusions  con- 
cerning them  are  somewhat  tentative,  but  are  probably  reasonable. 

The  expressions  used  for  and  the  fairing  seem  to  have  only  a 

minor  effect  on  calculated  velocity  profiles.  The  value  of  6*  can 
have  more  of  an  effect,  and  its  streamwise  variation  is  of  signifi- 
cance. If  6*  increases  too  slowly  in  the  streamwise  direction,  the 
velocity  gradient  near  the  wall  tends  to  remain  too  high.  Since  the 
value  of  is  of  some  significance,  a fairing  that  allows  the  eddy 
viscosity  to  remain  essentially  equal  to  for  a short  distance 
seems  appropriate.  Otherwise,  the  value  of  €+  may  be  effectively 
lost  since  the  calculations  may  fail  to  distinguish  where  the 
expression  stops  and  the  fairing  begins. 

Eq.  (5-31)  was  Included  in  the  formulation  at  one  time,  but  ten- 
ded to  reinforce  an  oscillation  that  was  occurring  in  the  solution, 
and  was  therefore  dropped.  The  oscillation  was  later  eliminated, 
but  use  of  Eq.  (5-31)  was  not  reconsidered.  For  the  fairing  used, 
the  value  of  £ * at  Y - \ depends  on  and  Kc  6 ^ • The  value 

of  6+is  increasing  very  rapidly  with  y for  this  portion  of  the 
profile,  and  the  value  at  in  reasonable  agreement  with 

Eq.  (5-32). 

Turbulent  wall  jets  typically  have  a negative  shear  stress  at 
/ — y ^ . The  explanation  for  this  is  that  the  turbulence  from 

the  outer  region  can  spill  over  into  the  inner  region,  and  vice 
versa,  but  the  large  eddies  in  the  outer  region  tend  to  dominate  the 


flow  at  y - y ' . Gartshore  and  Newman  (Ref.  19)  use  the  empirical 

expression 


for  the  magnitude  of  2”  at  Y/km  • Kind  (Ref.  18)  uses  the  same 

expression  but  with  a modification  for  curvature.  According  to  the 
above  expression,  is  of  the  order  of  20%  of  the  maximum  negative 
shear  stress  in  the  outer  region,  and  thus  can  have  a magnitude  com- 
parable to  or  greater  than  the  wall  shear  stress.  Gartshore  and 
Newman  found  that  using  =.  o in  their  Integral  method  resulted 

in  y being  too  small,  i.e.,  the  velocity  maximum  remaining  too 
close  to  the  wall,  and  a similar  result  was  obtained  in  the  present 
method  using  a conventional  eddy  viscosity  model.  A method  to  allow 
to  be  negative  was  therefore  devised. 

Following  the  usual  mixing  length  analysis  (Ref.  29),  but  re- 
taining second  derivatives  in  the  Taylor  series  and  applying  the 
analysis  near  the  velocity  maximum  gives 


where  is  the  mixing  length.  By  analogy  with  usual  mixing  length 
procedure,  the  shear  stress  is  taken  as 


72 


This  result  only  applies  near  the  velocity  maximum  and  it  is  not  de- 
sirable to  retain  the  second  derivative,  so  the  approximate  velocity 
profiles  given  by  Eq.  (5-15)  are  used  to  evaluate  the  term.  From 
Eq.  (5— 15b)  assuming  a profile  without  a minimum  velocity 
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dzv. 

_ Zb 
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A-O 

where  2b  *«2.1n2  * 1.387.  Using  the  mixing  length  expression 

= J'ay’J  } (5-16)  for  VC*  , the  above  expression  for 

evaluated  at  y=  L0  » and  <T  — ,*/LQ  as  typical  values  gives 
. %l_0  . Therefore,  the  shear  stress  correction  term  is 
taken  as 


£ 

1° 


- -V  (/+  £ +) 


. /3?7 


[/-/.3ff7(^=/]  (a-a.) 


This  yields  2^  - , which  is  in  reasonable  agreement  with 
experiment.  The  maximum  negative  shear  stress  is  unchanged  by  the 
above  term  since  the  term  vanishes  at  about  X^v>+  , $5"  L0' 

Use  of  Eq.  (5-15a)  gives 


dU-  _ A JtL 

*y  Y,~  \ >W  ~ y^  C y ) 

2l!YL  - ^ /_y  \H-i  kj(w-/)u  f y^N1 

yz  \Y^J  ~ y*-  \ yJ 
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Neither  of  these  expressions  match  the  outer  expressions  at  y~  ^ 


Also,  the  form  is  different,  and 


DHl 


does  not  become  zero  at  a 


convenient  point  as  it  does  for  the  outer  portion.  The  result  for 


1 yx 


is  therefore  used  as  a guide,  but  is  altered  so  that  it 


matches  the  outer  expression  to  give 


-pr  = -v 


'38 7 , y^  ,2.  . 

4*  ( y)  (&~  Me) 


the  ^ — j term  is  replaced  by  a linear  fairing  from  the  wall  to  the 
somewhat  arbitrarily  selected  point  / = , 6 y . Reasonable  variations 
in  this  matching  point  have  little  effect  on  calculated  results.  Re- 
sults of  the  above  analysis  are  added  to  the  conventional  eddy  vis- 
cosity model,  giving  Eq.  (5-14),  where 


X*  = 

r = 

X*  = 


• /3S7 


./387 

7/o 


7/, 


(5-34a) 

(5-34b) 

■T1 

(5-34c) 

(5-34d) 


Eq.  (5-14)  along  with  Eq.  (5-34)  is  used  for  profiles  with  or  without 
a minimum  velocity.  No  direct  modification  is  made  for  curvature. 

The  equations  are  solved  by  a finite  difference  method  based  on 
the  Keller,  Cebeci  method  (Ref.  14).  This  method  applies  to  a system 
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of  first  order  differential  equations,  and  the  present  equations  are 
rewritten  as  such  a system  by  defining 


G - 

(5-35) 
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//  = 

D<S 

(5-36) 
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With  these  definitions,  Eqs.  (5-8,  9,  10,  11,  14)  become 
Momentum  x 


y \_c2  0^]  - 4 &J  af  O^F) 


ZL £ 


[h-  -fc- 
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Momentum  y 


ut 6 1 
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Inner  Boundary  Conditions 


r = o t (3-0  at  ^ = o 


Outer  Boundary  Conditions 


(5-37) 


(5-38) 


(5-39) 


G - / , Cz  - C y - 7> 


(5-40) 
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The  finite  difference  procedure  applies  to  the  arbitrary  grid 
" JSLot  ) + A1  ~ Zt  ..f  M 


7,  = 1-, + 1, 


r- 


'/ 


/=  J 


Values  at  a grid  point  t't  ‘V .)  are  denoted  as  F~f* , and  the  nota- 

'■>1  V / 


tlon 


- i(t,+  L-,)  , V^.K  ~ i (%■  *■  %, ) 

is  used  for  averages  of  grid  points.  Averages  are  taken  as 
I f r f-M-t  \ ^ J_  / /—  M 


The  finite  difference  procedure  uses  first  order  centered  dif- 
ference formulas  for  derivatives,  and  average  values  for  algebraic 
terms.  Such  a procedure  is  Inherently  stable  in  the  sense  that  no 
extraneous  solutions  to  the  difference  equations  can  be  present  and 
grow  as  calculations  proceed  downstream.  Other  types  of  instabilit- 
ies are  possible,  but  can  be  avoided.  Streamwise  derivatives  are 
taken  as 


( wf  - ) /t 


An  equivalent  expression  for  normal  derivatives  is 
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9M / 


c V*-  *£*)/% 


Because  of  the  form  of  the  equations,  some  simplification  without 
loss  of  accuracy  results  by  taking 


37 


/t- 


Eqs.  (5-35,  36)  are  taken  as 


(rf-KW-  ■ 

(«;-  s, :.)//?■  - 


This  is  the  form  used  by  Keller,  Cebeci.  It  results  in  no  loss  of 
accuracy  since  satisfying  the  above  equations  at  each  ^ also  re- 
sults in  averages  at  ^ being  satisfied.  Eq.  (5-38)  is  taken  as 

< %■  - > /4-  = 'i ( 

This  equation  could  also  be  written  for  ^ Instead  of  , 

and  initially  it  was.  The  difference  in  the  solution  was  insignifi- 
cant, so  the  simpler  form  above  is  used.  Eq.  (5-37)  is  taken  as 
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In  the  above  form,  the  shear  stress  term  is  the  only  '/  derivative, 
and  is  treated  in  the  form 


Because  of  the  other  terms  in  the  equation,  it  would  seem  better  to 
treat  it  in  the  form 


as  Keller,  Cebeci  do.  However,  the  second  form  results  in  a dif- 
ficulty if  the  calculations  begin  at  the  slot,  as  they  do  in  the 
present  solution.  As  discussed  later,  the  initial  wall  jet  profile 
is  taken  as  a fully  developed  slot  profile  plus  the  known  upstream 
boundary  layer  profile.  For  fully  developed  slot  flow,  the  shear 
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stress  decreases  linearly  from  the  wall  and  Its  slope  depends  on  the 
pressure  gradient  in  the  slot.  The  wall  jet  adapts  to  the  external 
pressure  gradient,  which  in  general  is  much  different  than  that  in 
the  slot,  so  the  shear  stress  gradient  near  the  wall  is  considerably 
different  for  the  first  downstream  point.  Use  of  the  second  finite 
difference  form  therefore  results  in  drastic  oscillations  as  the  cal- 
culations try  to  make  the  average  shear  stress  gradient  correct. 

The  first  form  avoids  this  difficulty.  It  would  be  possible  to 
switch  from  the  first  to  the  second  form  at  some  point  downstream  of 
the  slot.  At  one  time  this  was  done,  but  it  resulted  in  an  insigni- 
ficant difference,  so  the  first  form  is  retained  throughout  for 
simplicity. 

In  the  calculations,  a solution  is  known  at  a station  f , 

J/A-  } 

and  the  equations  are  used  to  find  the  solution  at  . Use  of  the 
finite  difference  expressions  results  in  the  above  system  of  non- 
linear algebraic  equations.  An  iterative  procedure  is  required  to 

solve  these  equations,  and  the  Keller,  Cebeci  procedure  is  used.  The 

• — (l) 

superscript  /y\  is  dropped,  and  the  l ilx  iteration  is  denoted  as  j~.  . 

__  ( t'-hi)  (i)  ( z) 

Also,  f~ . ■ p , -h£F,  with  similar  expressions  for  the  other 

1 1 °r7 

variables.  The  is  a perturbation  term  which  is  solved  for  in 

r 

each  iteration.  Substituting  the  above  expressions  into  the  finite 
difference  equations  and  neglecting  products  of  perturbation  terms. 
From  Eq.  (5-35) 

<Fci>  P cr cP‘  * &L  = 9?' 

*7  bGr  ~Shr<  + -f  h 


(5-41) 
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where  <2 , — — p,  i~  P • + 

0r'tL  1 


t;  41 


From  Eq.  (5-36) 
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(5-42) 
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From  Eq.  (5-38) 
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From  Eq.  (5-37) 
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Eq.  (5-44)  neglects  terms  resulting  from  &£*  since  there  is  no 
simple  way  to  handle  such  terms.  Keller,  Cebeci  also  neglected  such 
terms  and  found  that  doing  so  did  not  damage  convergence. 

The  boundary  conditions  (5-39)  and  (5-40)  are  enforced  for  each 
iteration.  Since  ^ •=  O , Eq.  (5-39)  results  in 

S O (5-45a) 

bG,  - O (5-45b) 


The  thickness  of  the  jet  at  a station  ^ is  determined  during  the 
solution.  A basic  grid  is  selcted,  and  the  outer  point  7j-  in  the 
grid  should  be  far  enough  from  the  wall  to  allow  for  the  maximum  jet 
thickness  that  will  be  encountered.  Over  much  of  the  jet  region,  the 
jet  will  be  much  thinner  than  this,  so  the  entire  grid  is  not  used  in 
all  calculations.  The  outer  boundary  conditions  are  handled  as 
follows.  A grid  point  6=  J~  , corresponding  approximately  to 

=z  /.2~  7^  ( , is  selected  as  the  outer  point.  Eq.  (5-12) 

or  (5-13)  is  applied  to  obtain  U*  and  at  ^ je  • Eq*  (5-40)  is 

enforced  at  f°r  each  iteration,  and  this  gives 

c>  — & (5— 46a) 


(5-46b) 


The  actual  thickness  ^ of  the  jet  may  be  less  than 
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The  system  of  equations  is  solved  by  a matrix  factorization  pro- 
cedure analogous  to  that  used  by  Keller,  Cebeci.  The  equations  are 
written  in  a particular  order  and  then  considered  in  groups  of  four 
to  obtain  the  desired  form.  The  first  group  consists  of  Eqs.  (5-45a) 
(5-45b),  (5-44)  for  j « 2,  and  (5-43)  for  j ■ 2,  in  that  order. 

The  second  group  (in  order)  consists  of  Eqs.  (5-41)  for  j * 2, 

(5-42)  for  j * 2,  (5-44)  for  j ■ 3,  and  (5-43)  for  j ■ 3. 
Subsequent  groups,  except  for  the  last  group,  are  analogous  to  the 
second  group,  but  with  j increased  by  one  for  each  group.  The 
last  group  (in  order)  consists  of  Eqs.  (5-41)  for  j = Jg,  (5-42) 
for  j ■ J , (5-46a) , and  (5-46b) . This  system  of  equations  in 
matrix  form  is  then 

[s"']  = [<?"’] 

where 
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The  system  O'"  11 V'j  •[■?•]  can  be  solved  efficiently 

by  matrix  factorization,  and  this  is  the  procedure  used  by  Keller, 
Cebecl.  Let 
where 
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where 


_ <e  i r\ 

T),  = H, 


to 


ti)  «■  I -0<‘>  r-<‘> 

^ = \ ~ B/-  ^ 


e/>=  c? 


The  (iXty)  matrices  } By  } Qf  * 

viously,  and  the  (tyx.*/)  matrices  X)  * 


(t) 


have  been  given  pre- 
(O 


and 


can  be 


calculated  sequentially  from  them  as  shown  above.  The  matrix  [ ~Dj  J 


C*) 

is  the  inverse  matrix  of  By  . The  matrix  inversion  is  performed 
using  a standard  matrix  factorization  technique  (Ref.  30).  It  is 
necessary  that  f\^  be  nonsingular,  and  the  equations  were  arrang- 
ed in  the  order  used  to  accomplish  this. 

With  0“'J]  a»d  O'*  J known,  the  matrix  equation  can 
be  written  as 

[X-'J  [£->]-  [?*"] 

and  solved  for  2 ^ J where  £ 2 ^ ^ 


With 


where 


the  solution  is 
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[ z"‘  ] known,  Ort]  is  found  from 
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The  solution  is 
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r <-(<)"]  r -M  _<■<■)  .w  (c\  n 

With  |_  o J known,  new  values  yly  } Qy  , Cfy  ) 

are  determined,  and  the  procedure  is  repeated  until  convergence  oc- 
curs. When  convergence  occurs,  the  solution  for  ^ is  known,  and 
calculations  proceed  to  j A(+ 1 . Convergence  criteria  are  taken  as 

/ i* ) £ * 

b A/|  and  £ becoming  less  than  some  prescribed  tolerance. 

This  is  equivalent  to  the  velocity  gradient  at  the  wall  and  the  in- 
tegral of  the  velocity  profile  squared  not  changing  significantly 
for  the  t -tk  iteration. 

Separation  is  determined  by  the  calculated  value  of  be- 
coming zero  or  negative.  For  the  initial  guess  for  the  velocity 
profile  at  a station  ^ /-/  ^ (i.e.,  /-f {<r)  is  positive.  Subse- 
quent iterations  give  — // ^ until  convergence 

. (iti)  lit1*1'* 

occurs  or  /■+  40*  If  n,  ^ O , iterations  terminate  and  sep- 

aration is  assumed  to  have  occurred. 

Wall  jet  calculations  begin  at  the  blowing  slot.  A velocity 
profile  for  the  flow  in  the  slot  is  assumed,  and  the  boundary  layer 
profile  above  the  slot  is  known  from  the  upstream  boundary  layer 
calculations.  The  flow  in  the  slot  is  assumed  to  be  a fully  devel- 
oped incompressible  channel  flow.  For  such  a flow  (Ref.  31) 


ut  _ s!Z 

^7  - 7*  ' constant 

Assuming  symmetry  about  the  slot  centerline,  this  gives 

= 'tux  C 1 - if) 


(5-47) 
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where  (t  * slot  thickness.  For  laminar  flow, 

and  the  velocity  profile  corresponding  to  Eq.  (5-47)  is 

'fc’io-  / y 


«=  ( I - ir) 


-v  /° 

If  a suitable  eddy  viscosity  distribution  is  known  for  turbulent 
flow,  a similar  procedure  can  be  used  to  obtain 
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The  eddy  viscosity  model  developed  by  Mei  and  Squire  (Refs.  32,  33) 
gives 
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This  expression  along  with  Eq.  (5-47)  gives 
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Integrating  Eq.  (5-48)  gives  G and  F.  All  terras  in  Eq.  (5-48)  are 
known  except  L/*  , which  must  be  prescribed.  Prescribing  thus 
prescribes  the  blowing  slot  values.  Eq.  (5-48)  can  be  integrated  in 
closed  form  (Refs.  32,  33)  but  for  the  present  purpose  it  is  more 
convenient  to  do  the  integration  numerically.  Either  a turbulent 
slot  flow  or  a laminar  slot  flow  can  be  prescribed.  For  the  laminar 
case,  the  eddy  viscosity  is  set  equal  to  zero.  For  all  calculations 
performed  so  far,  a turbulent  slot  flow  has  been  assumed. 

Blowing  slot  values  are  normally  given  in  terms  of  a momentum 
coefficient  and  a mass  flow  coefficient  CyA  , where 

(yi  - -f  ^ f d-y  = «-«> 

X C o 
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p.  j[  udy 
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(5-50) 


In  Eqs.  (5-49,50)  is  the  density  of  the  flow  in  the  slot.  Since 
the  wall  jet  analysis  is  for  incompressible  flow,  it  is  assumed  that 
f*  • Eq.  (5-49)  is  the  appropriate  form  for  the  case  where  the 
velocity  in  the  slot  is  not  uniform.  For  most  cases,  however,  ex- 
perimental values  of  are  based  on  an  average  velocity.  An  av- 
erage velocity  &%)/£  for  the  slot  is 

c 
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Defining  a momentum  coefficient  CL,.  , in  terms  of  the  average 

Av£ 


velocity  gives 
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(5-51) 


The  values  predicted  by  Eqs.  (5-49)  and  (5-51)  are  reasonably  close 

but  not  identical.  For  Kind's  Flow  II,  assuming  a turbulent  slot 

flow,  the  value  of  C ,,  is  about  102  higher  than  £. ,, 

s*-  Avc 

With  (A^  for  the  slot  prescribed,  the  velocity  profile  and  thus 
and  are  known.  The  upstream  boundary  layer  profile  is 
also  known,  but  in  terms  of  a different  coordinate  system.  Letting 
the  subscript  ( ) denote  quantities  in  the  upstream  boundary 
layer  coordinates  (evaluated  at  the  slot) , the  values  in  the  wall 
jet  coordinates  become 
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Neither  the  upstream  boundary  layer  calculations  nor  the  slot  flow 
analysis  include  normal  pressure  gradients,  so  the  initial  pressure 
distribution  is  taken  as  C ^ = / — Ct*  • For  the  form  used  for 
the  finite  difference  expressions,  this  does  not  cause  any  difficulty, 
and  neither  does  the  fact  that  /-/  is  discontinuous  at  ^ =•  / . 

Large  gradients  are  associated  with  the  sharp  velocity  minimum 
at  the  slot,  so  a finer  grid  is  necessary  near  the  slot.  From  re- 
sults presented  by  Keller,  Cebeci  (Ref.  14),  40  ^ grid  points  seemed 
a reasonable  number  for  a conventional  boundary  layer  profile.  The 
shape  of  the  initial  wall  jet  profile  requires  tripling  this  number, 
so  a maximum  of  121  grid  points  was  selected.  The  //  grid  for 
the  initial  calculations  is  taken  as  follows.  For  the  region 
the  Cebeci,  Smith  grid  =-  is  used,  but  the  value  of  K 

is  altered  occasionally.  For  ^ / , the  grid  is  a mirror  im- 
age of  the  grid  for  Selected  points  in  the  upstream  boundary 

layer  grid  are  used  for  the  remaining  wall  jet  grid  points.  The 
first  streamwlse  step  size  is  taken  as  1/3  of  the  slot  thickness, 
and  the  second  step  size  is  taken  as  1 slot  thickness.  After  the 
solution  is  obtained  at  this  location,  the  grid  is  altered  to  a more 
suitable  one  for  downstream  calculations.  The  Cebeci,  Smith  grid 

P . — % • is  used,  again  with  £ altered  occasionally.  The 

• ■ j / " / 

streamwlse  step  size  can  be  varied  as  desired,  but  unless  output  at 
specific  points  is  desired  the  step  size  is  taken  as  follows.  A max- 
imum step  size  (such  as  5 slot  thicknesses)  is  prescribed.  This  is 

taken  as  the  step  size  unless  it  results  in  changing  by  more  than 

£ 

.1  between  successive  locations.  If  Uq  changes  by  more  than  .1,  the 
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step  size  is  reduced  so  as  to  make  the  change  in  approximately 
.1.  However,  the  step  size  is  not  reduced  to  less  than  half  the  pre- 
scribed maximum  value.  The  values  of  K were  selected  on  a trial 
and  error  basis,  and  general  recommendations  on  how  to  pick  values 
would  require  additional  study.  Calculations  were  performed  with 
different  grid  parameters  to  insure  that  the  grid  used  is  adequate 
for  the  present  calculations. 

The  grid  used  for  the  wall  jet  calculations  is  quite  coarse 
compared  to  that  used  for  the  boundary  layer  calculations.  Prelim- 
inary calculations  indicate  that  an  even  coarser  grid  would  be  ac- 
ceptable, although  oscillations  in  the  solution  very  near  the  wall 
can  result  if  the  grid  is  too  coarse  near  the  wall.  Convergence 
of  the  finite  difference  method  was  found  to  be  fast  and  reliable. 
However,  a convergence  anomaly  existed  during  much  of  the  develop- 
ment, and  its  presence  was  not  directly  obvious.  The  initial  guess 
for  the  velocity  profile  at  a new  downstream  location  is  taken  as 
the  solution  from  the  previous  location,  and  this  seems  to  be  satis- 
factory. However,  unless  some  precaution  is  used,  this  can  result  in 
the  initial  guess  for  — from  Eq.  (5-26)  — being  too  small  while 
— from  Eq.  (5-25)  — increases  too  rapidly.  This  results  in 
a qualitatively  poor  initial  eddy  viscosity  distribution,  and  subse- 
quent iterations  do  not  seem  to  be  able  to  compensate  even  though 
convergence  appears  to  be  normal.  Providing  a qualitatively  reason- 
able Initial  eddy  viscosity  distribution  eliminates  the  problem,  and 
the  final  solution  is  not  sensitive  to  the  initial  guess  used. 
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Summary  Of  Wall  Jet  Calculations 

Wall  jet  calculations  begin  at  the  blowing  slot.  A value  of 
Ut  for  the  slot  is  assumed , and  the  velocity  profile  for  the  slot 
is  determined  from  Eq.  (5-48).  The  rest  of  the  Initial  velocity  pro- 
file is  determined  from  the  known  upstream  boundary  layer  profile  and 

is  given  by  Eq.  (5-52).  The  slot  location  is  denoted  as  ? . The 

Ji 

solution  proceeds  downstream  from  the  slot.  The  solution  at  £ and 

is  determined  using  a finite  difference  grid  appropriate  to  a 
velocity  profile  with  a sharp  minimum.  After  the  solution  at  ^ is 
obtained,  the  grid  is  altered  to  a form  more  appropriate  for  down- 
stream calculations,  and  the  solution  at  the  new  grid  points  is  ob- 
tained by  interpolation.  The  solution  then  proceeds  downstream.  The 
initial  wall  jet  calculations  are  for  a profile  with  a relative  min- 
imum velocity.  For  such  a profile,  the  outer  edge  conditions  (at  a 
point  corresponding  to  the  outer  edge  of  the  upstream  boundary  layer) 
are  determined  from  Eq.  (5-13).  To  compensate  for  the  lack  of  a cur- 
vature correction  for  Eq.  (5-13),  the  normal  pressure  gradient  in  the 
wall  jet  equations  is  set  equal  to  zero  above  the  velocity  minimum. 
The  eddy  viscosity  distribution  is  determined  from  Eqs.  (5-25,  26, 

33,  21,  22,  23). 

Calculations  retain  the  minimum  velocity  until  they  can  switch 
smoothly  to  a profile  without  a minimum.  Far  profiles  without  a vel- 
ocity minimum,  the  outer  edge  conditions  are  determined  from  Eq. 
(5-12).  The  eddy  viscosity  distribution  is  determined  from  Eqs. 
(5-25,  26,  33,  19,  20).  Calculations  proceed  downstream  until  sep- 
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aration  is  determined  by  the  wall  shear  stress  becoming  zero  or  neg- 
ative. 

For  all  wall  jet  calculations,  the  outer  boundary  conditions 
are  applied  at  a grid  point  corresponding  approximately  to  /*Z  ’vj 
The  actual  jet  thickness  is  determined  during  the  solution  and  may 

i 

be  less  than  /.  Z.  The  solution  for  'j  ? Tg  is  taken  as 

F - F7_  +-  (/t]-^}je)  j <3-/  , H=0)  Ci=^Jhe  initial  8uess  for  //> 

at  a station  ^ is  taken  as  the  corresponding  values  at  • 

The  Initial  guess  for  C P is  taken  as  Cp  + -6  Cp  where  4 C,p 

Al  J~M-i 

is  determined  so  as  to  give  the  correct  outer  edge  pressure.  After 
the  first  iteration  (only)  at  a station  , the  old  and  new  solu- 
tions are  averaged.  This  did  not  seem  to  be  necessary,  but  helped 
reduce  over-shooting  of  the  next  iteration.  A reasonable  eddy  vis- 
cosity distribution  is  provided  for  the  first  iteration. 

An  iterative  procedure  is  required  to  determine  the  correct  blow- 
ing slot  values.  Such  a procedure  should  not  be  difficult  to  set  up, 
but  at  present  calculations  are  simply  performed  for  a given  number 


of  assumed  slot  conditions. 
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6.  RESULTS 

Calculations  are  performed  for  Kind's  Flow  II  (Refs.  4,  5,  15). 
Briefly,  Kind's  airfoil  is  a 20%  thick  ellipse  with  the  trailing 
edge  replaced  by  a circular  section  with  a radius  of  9/16  inches. 

The  chord  length  is  14  5/8  inches,  and  the  blowing  slot  is  located 
at  the  location  where  the  modified  trailing  edge  joins  the  airfoil. 
Kind  performed  downwash  corrections  by  comparing  calculated  and  mea- 
sured pressure  distributions  and  estimates  the  accuracy  of  the  re- 
sulting effective  two  dimensional  angle  of  attack  as  + 1/2°.  Kind's 
experimental  values  for  Flow  II  are  ^<2.  on  ■ 750000,  * -0.7°, 

Cjl  - 1.82,  .055,  <£*«=  .0012,  £*=  .0378. 

Using  this  input,  calculations  are  performed  as  described  pre- 
viously. Results  of  the  wall  jet  calculations  are  given  for  two 
different  sets  of  blowing  slot  conditions  in  order  to  illustrate  the 
effects  of  the  assumed  slot  conditions.  Results  of  the  present  cal- 
culations are  shown  in  Figs.  6 through  15  and  are  discussed  below. 

Fig.  6 shows  the  calculated  potential  flow  surface  pressure  dis- 
tribution. Unfortunately,  Kind  does  not  give  the  measured  distribu- 
tion for  Flow  II,  so  a comparison  is  impossible.  The  same  potential 
flow  program  used  here  has  been  used  in  other  applications,  such  as 
in  Ref.  3,  where  calculated  values  were  found  to  be  in  good  agreement 
with  experiment  for  a cambered  ellipse  with  a modified  trailing  edge. 

Fig.  7 shows  calculated  and  measured  boundary  layer  profiles  up- 
stream of  the  blowing  slot,  and  the  agreement  is  good.  The  measured 
profile  is  .055C  upstream  of  the  slot,  while  the  calculated  profile 
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is  .007C  upstream  and  is  the  profile  used  for  the  wall  jet  calcula- 
tions. The  difference  in  location  could  cause  some  difference  in 
the  profiles,  but  there  seems  to  be  little  upstream  influence  of  the 
blowing  slot  for  this  flow  (Potential  flow  calculations  give  a 
at  the  slot  of  about  -2.3,  while  Kind's  measured  values  give  about 
-2.5  at  the  slot  exit).  Other  potential  sources  of  difference  are 
the  relatively  low  Reynolds  number  and  the  fact  that  Kind  used  two 
trip  wires  on  both  the  upper  and  lower  surfaces  of  the  airfoil.  As 
mentioned  previously,  the  boundary  layer  eddy  viscosity  expressions 
were  not  updated  by  the  newer  expressions  for  low  Reynolds  number. 

In  calculations,  transition  was  assumed  at  the  first  trip  wire,  and 
the  second  wire  was  ignored. 

Fig.  8 shows  calculated  and  measured  velocities  at  the  outer 
edge  of  the  wall  jet,  along  with  the  calculated  potential  flow 
velocity.  Calculated  values  shown  are  for  S*-  " .054,  but  values 
for  9*-  .046  are  virtually  the  same.  Calculations  include  a vel- 
ocity minimum  up  to  ~ slot**  *011,  and  the  dashed  line  indicates 
the  calculated  value  of  • Calculated  values  are  In  excellent 

agreement  with  experiment  up  to  about  j-  % slot  * .05,  and  then  de- 
crease in  agreement,  but  are  still  useable.  The  discrepancy  begins 
at  about  the  same  place  where  the  normal  velocity  becomes  comparable 
to  the  streamwise  velocity,  and  other  difficulties  in  the  wall  jet 
calculations  begin  at  about  the  same  location.  Eqs.  (5-1)  and  (5-2) 
give  virtually  the  same  result  up  to  | 5 Lor-  m .05,  with  Eq. 

(5-2)  giving  some  improvement  downstream.  The  discrepancy  in  the 
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outer  edge  velocity  forces  some  discrepancy  in  calculated  velocity 
profiles  (G)  also.  Fortunately,  however,  it  seems  that  the  effects 
on  the  velocity  ( G>  ) in  the  inner  portion  of  the  jet  are  not 
great.  This  is  discussed  later. 

Fig.  9 shows  calculated  and  measured  pressure  distributions  in 
the  wall  jet  region.  The  experimental  outer  edge  pressures  were 
calculated  by  Kind  using  his  measured  velocities  and  the  assumption 
that  the  total  pressure  remained  constant  along  a streamline.  Cal- 
culated values  are  shown  for  = .054.  Values  for  ” .046 

are  slightly  different  but  are  close  enough  that  the  curve  is  some- 
what difficult  to  read  if  they  are  included.  Calculated  outer  edge 
pressures  are  for  the  point  in  the  wall  jet  finite  difference  grid 
where  the  outer  boundary  conditions  are  applied. 

Calculated  and  measured  pressures  are  in  reasonable  agreement 
with  each  other.  The  calculated  wall  pressure  gradient  is  in  good 
agreement  with  experiment  for  portions  of  the  jet  region,  but  is 
higher  in  places  and  lower  in  others.  The  difference  is  of  the  order 
of  30Z  for  some  of  the  latter  portion  of  the  jet,  and  this  has  some 
observable  effects  on  the  velocity  profiles.  However,  the  results 
are  still  reasonably  good.  In  contrast  to  this  is  the  case  indicated 
for  f ( f ) ”0.  The  calculations  for  this  case  are  identical  ex- 

cept that  the  approximate  streamline  curvature  correction  is  not  in- 
cluded. Up  to  about  % - $Acrm  *05,  the  results  for  this  case  are 
indistinguishable  from  the  previous  case.  However,  the  results  then 
begin  to  diverge  at  an  Increasing  rate.  As  the  wall  pressure  gra- 
dient remains  too  mild,  the  maximum  velocity  in  the  wall  jet  remains 
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too  large,  and  these  effects  reinforce  each  other.  The  end  result  of 
this  is  that  the  calculations  become  useless. 

The  rapid  change  in  wall  pressure  near  the  slot  is  a result  of 

fairing  in  the  radius  of  curvature  over  a short  distance.  A some- 

what longer  fairing  could  also  be  used,  but  should  be  completed  by 
about  £ — .01  to  match  the  measured  wall  pressure  values  in 

.01  < £ — 3st-or  < ° * Such  a fairing  is  still  quite 

rapid  to  apply  to  the  outer  boundary  conditions  at  the  outer  edge 
of  the  upstream  boundary  layer,  so  the  procedure  described  in  Section 
5 was  adopted.  The  actual  fairing  used  does  not  seem  significant  in 
the  downstream  calculations. 

The  pressure  in  the  separated  region  is  in  reasonable  agree- 
ment. Kind's  measured  value  for  C. 2 ^ is  .65,  the  calculated  value 
for  the  lower  surface  is  .47,  and  the  calculated  value  for  the  upper 
surface  is  .27  for  “ .054  ( — Z for  m .046).  The  calcul- 

ated value  is  taken  as  the  value  obtained  on  the  last  iteration, 
l.e.,  the  iteration  that  determined  a separated  profile.  However, 
it  is  somewhat  difficult  to  get  a precise  value  for  the  wall  jet  sep- 
aration pressure.  The  pressure  gradient  is  very  large,  and  a small 
change  in  the  location  of  the  separation  point  can  make  a significant 
change  in  the  calculated  pressure.  This  is  complicated  by  the  ex- 
pression used  for  -f  ( j ) being  crude,  and  can  also  be  complicated 
by  using  a fairly  large  step  size  as  was  done  in  the  present  calcul- 
ations . 

Figs.  10  through  15  show  calculated  and  measured  wall  jet  pro- 
files. The  values  of  — J ^given  are  for  the  corresponding  sta- 
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tions  In  the  wall  jet  grid.  In  some  cases,  the  stations  do  not  cor- 
respond exactly  with  the  location  of  the  measured  profile,  but  the 
difference  is  small.  Figs.  10  through  12  are  for  .054 

( Cfr  ■ .0055)  and  Figs.  13  through  15  are  for  ■ .046 
( “ .0051).  In  Fig.  12,  the  calculated  and  measured  profiles  on 

the  right  are  at  somewhat  different  locations. 

The  agreement  with  experiment  is  good  for  either  set  of  profiles, 
especially  since  measured  profiles  and  pressures  are  not  used  in  the 
calculations.  Which  set  provides  the  best  agreement  is  somewhat  am- 
biguous since  one  set  is  better  in  some  ways  while  the  other  is 
better  in  others.  Comparing  Figs.  10  and  13  shows  that  the  solution 
for  ■ .054  gives  better  agreement  for  station  9,  while  the 
solution  for  *>  .046  gives  better  agreement  for  the  inner  region 

at  station  12.  In  either  case,  however,  the  agreement  between  mea- 
sured and  calculated  values  is  within  about  5Z.  For  either  set  of 
profiles,  the  agreement  at  stations  15  and  18  is  essentially  the  same 
as  the  agreement  at  station  12,  although  at  station  18  there  is  a 
slight  tendency  for  the  maximum  velocity  to  be  too  large.  Direct  com- 
parison with  experiment  of  either  profile  at  station  22  is  complicated 
by  the  fact  that  the  calculated  and  measured  outer  edge  velocities 
differ  by  about  10Z.  The  solid  lines  represent  the  calculated  pro- 
files, while  the  dashed  lines  represent  the  profiles  multiplied  by  the 
ratio  of  calculated  to  measured  outer  edge  velocity.  The  dashed  lines 
are  thus  probably  more  representative  of  the  velocity  for  the  inner 
region.  If  so,  the  solution  at  station  22  tends  to  follow  the  sol- 
ution at  station  18  except  that  the  velocity  at  station  22  remains  too 
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high  very  near  the  wall.  For  ■ .046,  separation  was  calculated  at 
station  25  ( ^SLt>Tm  .0681)  and  for  .054  separation  was  cal- 


culated at  station  26  ( J ^ ■ .0706).  Kind  experimentally  de- 


file given  in  Fig.  12  (on  the  right)  is  for  %SLbTm  .0706.  Thus 
the  calculated  and  measured  separation  points  agree  well,  but  there 
is  no  calculated  profile  to  compare  to  experiment.  The  calculated 
profile  in  Fig.  12  is  thus  for  station  25  (the  last  profile  deter- 
mined), but  the  step  size  is  relatively  large  and  agreement  should 
not  be  expected.  However,  it  is  of  some  interest  to  indicate  how 
rapidly  the  jet  is  changing,  so  the  profiles  are  included.  The  cal- 
culated profile  is  multiplied  by  the  ratio  of  the  calculated  to  mea- 
sured velocity. 

Without  additional  comparisons  with  experiment,  conclusions 
drawn  from  the  above  observations  are  somewhat  tentative.  However, 
during  (numerous!)  development  calculations,  some  results  occurred 
with  such  monotonous  regularity  that  some  effects  seem  clear.  All 
calculations  were  started  with  a turbulent  slot  profile.  The  vel- 
ocity profile  at  station  9 is  only  sensitive  to  the  assumed  CAL  and 
to  the  value  of  [Eq.  (5-19)  or  (5-21)).  Variations  in  the 

radius  of  curvature  fairing  and  inner  eddy  viscosity  model  have  little 
effect  on  the  portion  of  the  profile  far  enough  from  the  wall  to  be 
compared  to  experiment.  Until  pressure  gradient  effects  begin  to  be- 
come severe  (at  about  station  18),  the  solution  depends  largely  on 
the  value  of  Kc  . Between  stations  9 and  12,  there  is  a strong 

tendency  for  the  maximum  velocity  to  decrease  too  slowly  and  the  thick- 


termined  separation  at 
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ness  of  the  outer  region  to  remain  too  small.  Also  in  this  region, 
the  value  of  tends  to  remain  nearly  constant.  The  curvature 
correction  accounts  for  most  of  the  change  in  in  this 

region,  and  before  it  was  incorporated  the  discrepancy  between  cal- 
culated and  measured  values  was  about  double  that  in  the  present  re- 
sults. This  may  be  an  indication  of  some  difficulty  in  the  expres- 
sion for  C ^ as  the  flow  passes  from  a favorable  to  an  adverse  pres- 
sure gradient.  From  station  12  to  station  18,  calculated  velocity 
profiles,  pressure  gradients  and  outer  edge  velocities  follow  mea- 
sured values  very  well.  Pressure  gradient  effects  begin  to  become 
severe  at  about  station  18,  and  tend  to  dominate  calculations  down- 
stream. For  the  calculations  for  ^ " .054,  for  example,  the  shear 
stress  variation  over  O £.  y ± .Z  is  about  5Z  for  station  15,  15Z 
for  station  18,  50Z  for  station  22,  300Z  for  station  24,  and  separa- 
tion occurs  at  station  26.  Thus  pressure  effects  are  comparable  to 
shear  stress  effects  for  a considerable  distance  upstream  of  separa- 
tion. If  the  wall  pressure  gradient  is  too  mild,  the  velocity  very 
near  the  wall  tends  (strongly)  to  remain  too  high.  Fig.  9 shows  that 
the  wall  pressure  gradient  is  somewhat  too  mild  between  stations  18 
and  22,  and  the  velocity  profiles  at  station  22  show  some  of  this 
effect.  Downstream  of  this  the  calculated  wall  pressure  gradient 
becomes  larger  than  the  measured  one,  and  this  tends  to  compensate 
to  some  degree.  Calculations  over  the  latter  portion  of  the  jet  re- 
quire a reasonable  wall  pressure  distribution  if  they  are  to  have  any 
hope  of  success. 


103 


7.  CONCLUSIONS 

A self-contained  analysis  for  arbitrary  circulation  controlled 
airfoils  in  incompressible  flow  is  developed  and  found  to  be  capable 
of  predictions  that  are  in  good  agreement  with  experiment.  The  an- 
alysis requires  calculations  for  potential  flow,  laminar  and  turbu- 
lent boundary  layer  flow,  and  turbulent  wall  jet  flow. 

The  turbulent  wall  jet  calculations  are  the  most  difficult  and 
are  subject  to  the  most  uncertainty.  Accurate  calculations  require 
reasonably  accurate  predictions  of  the  wall  pressure  distribution  in 
the  wall  jet  region.  Conventional  boundary  layer  approximations  pro- 
vide an  accurate  wall  pressure  distribution  for  much  of  the  wall  jet 
region,  but  require  modification  for  the  latter  portion  of  the  jet. 

An  exact  correction  to  account  for  this  seems  untenable  since  such  a 
correction  tends  to  produce  stability  problems  in  the  calculations. 
Such  stability  problems  might  well  be  expected  since  the  correction 
involves  second  order  streamwise  derivatives,  which  do  not  fit  into 
boundary  layer  theory.  A plausible  approximate  correction  based  on 
average  values  seems  adequate  and  workable.  Such  a correction  is 
used,  and  suggestions  for  improving  it  are  offered. 

The  shear  stress  model  for  the  wall  jet  is  subject  to  some  un- 
certainty because  of  a lack  of  detailed  experimental  data,  or  in  some 
cases  because  of  conflicting  experimental  data.  However,  the  model 
does  seem  reasonable  and  is  formulated  using  available  information. 

It  is  possible  that  additional  comparison  of  calculated  and  experi- 
mental results  would  lead  to  improvement  of  the  model.  Such  compar- 
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isons  would  probably  be  more  helpful  if  the  flow  were  subject  to  less 
severe  curvature  and  pressure  gradients  than  in  the  present  calcula- 
tions. 


i 
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APPENDIX:  COMPUTER  PROGRAM 

This  appendix  contains  a listing  of  the  computer  program  and 
a description  of  the  program  input  and  output.  The  program  consists 
of  the  potential  flow  program  described  in  Section  3,  the  boundary 
layer  program  described  in  Section  4,  and  the  wall  jet  program  de- 
scribed in  Section  5.  Only  two  dimensional  flows  ( € ■ 0)  without 
suction  ( (fur  “ 0)  are  considered  here,  and  the  boundary  layer  pro- 
gram is  set  up  for  this  type  of  flow. 

The  program  operates  as  follows.  With  the  input  given  (Table  Al), 
the  potential  flow  is  calculated  and  converted  to  the  form  required 
for  boundary  layer  and  wall  jet  use.  Boundary  layer  calculations  are 
performed  first  for  the  lower  surface  and  then  for  the  upper  surface. 
Wall  jet  calculations  are  then  performed  for  a prescribed  number  of 
assumed  slot  conditions.  The  initial  assumed  value  of  Uf£  for  the 
slot,  and  the  increment  ACi'tf  for  subsequent  assumed  values,  are 
input.  Drag  and  the  change  in  Cg  caused  by  the  pressure  variation 
over  the  wall  jet  region  could  be  calculated,  but  are  not  at  present. 

The  input  (Table  Al)  consists  of  three  subroutine  and  three  data 
cards.  The  subroutines  provide  the  airfoil  geometry  as  follows. 
Subroutine  BODY  provides  the  coordinates  ( y<z  ) of  the  airfoil  as 
required  for  the  Theodorsen  method.  Subroutine  SLOPE  provides  the 

c// 

values  of  — required  to  convert  to  curvilinear  coordinates.  For 

ofX  c . 

convenience  in  cases  where  analytic  expressions  for  — are  not 
available,  an  option  in  the  program  allows  values  of  <^)c/o^Xc to  be 
determined  numerically.  In  this  case,  subroutine  SLOPE  is  ignored. 


and  — is  automatically  calculated  numerically.  An  input  Integer 
tells  the  program  which  procedure  to  use.  Subroutine  RADIUS  provides 


the  surface  radius  of  curvature  In  the  wall  jet  region.  The  sub- 
routines for  Kind's  airfoil  are  Included  In  the  program  listing.  The 
three  data  cards  provide  the  rest  of  the  Input,  Including  specifying 
options  in  the  program.  The  form  required  for  the  data  cards  can  be 
determined  by  inspection  of  the  READ  statements  on  the  first  page  of 
the  program  listing. 

Output  from  the  program  consists  of  the  potential  flow  output 
(Table  A2),  the  boundary  layer  output  (Table  A3),  and  the  wall  jet 
output  (Table  A4) . The  numbering  system  is  different  in  the  three 
sets  of  output.  For  the  potential  flow  output,  points  are  numbered 
in  a counter-clockwise  direction,  with  point  1 being  the  trailing 
edge.  For  the  boundary  layer  output,  point  numbers  are  referenced 
from  the  front  stagnation  point.  Point  1 is  the  first  output  point 
downstream  of  the  stagnation  point  along  the  upper  surface  of  the 
airfoil.  The  remaining  upper  surface  points  are  numbered  consecu- 
tively until  the  rear  potential  flow  stagnation  point  is  reached. 

I 

If  the  last  upper  surface  point  is  denoted  by  jj.  , the  first  point 
downstream  of  the  front  stagnation  point  along  the  lower  surface  of 

t 

the  airfoil  is  denoted  by  + I ) • The  remaining  lower  surface 

points  are  numbered  consecutively.  The  wall  jet  points  are  numbered 
consecutively,  with  point  1 being  at  the  blowing  slot. 


125 


TABLE  Al:  COMPUTER  PROGRAM  INPUT 


Computer 

Symbol 

Quantity 
& Units 

Comments 

CL 

S’ 

Lift  coefficient 

ALPHA 

o<  <degrees) 

Angle  of  attack 

REC 

RGqo 

Free  stream  Reynolds  number 

NEPS 

Ne 

Number  of  output  points 

(must  be  even) . 

Recommended:  AO  — L.  80 

NBPTS 

nb 

Number  of  body  coordinates  to 

be  Input  (must  be  even) . 

Recommended : 200  N^  6.  400 

KUTTA 

Kt/ 

THjj  ■ 0 if  Input  is  to  be 

used. 

Ku  ■ 1 if  Kutta  condition  is 

to  be  used. 

NAEQN 

na 

N ■ 0 if  airfoil  slopes  are  to 
A 

be  calculated  numerically. 

N.  « 1 if  subroutine  SLOPE  is 
A 

to  be  used. 

AK 

K 

Turbulent  y grid  parameters 

HI 

hl 

i 

See  Table^JV  Ref.  8 , p>?  6? 

NT 

N 

l 

(continued  next  page) 
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CHORD 

c ,(£«et)  \ 

Chord 

BT 

XCJET 

b <feet) 

X 

C8lot 

See  Fig.  3 

Cartesian  coordinate  of 
bloving  slot 

SLOTT 

Dimensionless  slot  thickness 

UTSIG 

(A 

Initial  value  of  CCr  for  the 
blowing  slot 

UTS  IN 

A Ur 

Increment  of  U.T  for  the 
blowing  slot 

STEPS 

L/t* 

Maximum  £ step  size  for  wall 
jet  grid 

NWJIT 

Number  of  times  wall  jet 
calculations  are  to  be 
performed.  Assumed  Is 

increme  ted  by  ^ Ut  each  time. 

LAMSL 

LAMSL  ■ 0 for  a turbulent  slot 
profile. 

LAMSL  “ 1 for  a laminar  slot 
profile 

Subroutine  BODY 

Subroutine  BODY  provides  values 

for  XX(J),  YY(J)  (j-l,...,Nfi)  , 

c and  b.  XX  — * x and  YY  -V  y 

c c 

In  Fig.  3.  All  dimensions  are 

(continued  on  next  page) 
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In  feet.  Points  are  nuabered 
In  counterclockwise  order  around 
the  airfoil.  The  point  J ■ 1 
must  be  the  trailing  edge 
(y„  * 0,  x > 0)  and  J-1/2N+1 
must  be  the  leading  edge 
(yc  - 0,  xc<  0).  It  is 
desirable  to  have  nearly  equal 
arc  lengths  between  points. 

Subroutine  SLOPE 

Subroutine  SLOPE  provides  the 

slopes  of  the  upper  and  lower 

airfoil  surfaces  as  functions  of 

x . 
c 

/ otic-  \ 

YPU  “ \ cf/ c ' upper  surface. 

YPL  * ) lower  surface. 

The  function  statement 
X - . 5* (XUL-XLL ) *XA+ . 5* (XUL+XLL ) 
must  be  included. 

Subroutine  RADIUS 

Subroutine  RADIUS  provides  R* 
for  a given  X,  where  X*^-]?  ^ 
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TABLE  A2:  POTENTIAL  FLOW  OUTPUT 


Computer 

Symbol 

Quantity 
& Units 

Comments 

PHI 

V 

(degrees) 

Angle  <p  defined  In  Fig.  3. 

X 

X 

c 

(feet) 

fSee  Fig.  3 

Y 

yc 

(feet) 

l 

U 

Eq.  (3-14) 

CP 

cp 

Eq.  (3-16) 

CL 

CX 

Lift  coefficient 

ALPHA 

oC 

(degrees) 

Angle  of  attack 

XSTAG 

X . 

stag 

(feet) 

Coordinates  of  front 

YSTAG 

ystag 

(feet) 

stagnation  point. 
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TABLE  A3:  BOUNDARY  LAYER  OUTPUT 


Computer 

Symbol 

Quantity 
& Units 

Comments 

X 

x (feet) 

Fig.  1 

XI 

r 

Eq.  (3-17a) 

BETA 

Eq.  (3-17b) 

UE 

a* 

Eq.  (3-14) 

XC 

xc  (feet) 

Fig.  3 

YC 

yc  (feet) 

Fig.  3 

DELTA 

£ (Inches) 

Eq.  (4-23) 

DELST 

£*  (inches) 

Eq.  (4-21) 

THETA 

O (inches) 

Eq.  (4-22) 

CF 

Cf 

Eq.  (4-20) 

H 

H 

H - &*/e 

REC 

Re  oo 

REDS 

ReS* 

Eq.  (4-24) 

RET 

1 

Reo 

© 

Eq.  (4-25) 

ITERATION 

> 

Maximum  is  20 

DFW,  DDS 

Iteration  tolerances 

For  laminar  flow,  DDS  ■ 0. 

(continued  next  page) 
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NEDGE 

NEDGE  oust  be  < NT  (Table  Al), 
If  NEDGE  - NT,  the  ^ grid 
should  be  changed  to  allow  for 
a thicker  boundary  layer. 

FW" 

t" 

Values  of  f"  Indicate 

V 

w 

appropriate  values  of  h^ 
(Ref.  8). 
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TABLE  A4:  WALL  JET  OUTPUT 


Computer 

Symbol 

Quantity 
& Units 

Comments 

CMU 

Eq.  (5-49) 

CM 

Eq.  (5-50) 

UAVE 

Average  velocity  In  bloving  slot* 

UMAX 

U* 

An 

Maximum  velocity  In  bloving  slot. 

HW 

H 

w 

XI 

5 k 

Eq.  (5-7) 

DELXI 

UE 

u? 

Eq.  (5-12a)  or  Eq.  (5-13) 

CPW 

C 2 
xor* 

CPE 

Eq.  (5-12b)  or  Eq.  (5-13) 

CP  (POT) 

CP* 

o! 

1 

TAUW 

tt 

ZUatiur/  (£*%*«,) 

DELTA 

A * 

GMAX 

Gm 

GMIN 

Gmin 

Graln  ” G f°r  profiles  vithout 

a minimum  velocity 

ITERATION 

Maximum  is  10 

ETAMAX 

% 

ETAMIN 

^)min 

DHW,  DCPW 

Iteration  tolerances 

JDELTA 

"y  - at  / - JDELTA 

non 
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DIMENSION  ZP(401  ) ,Z0(401)  ,L(  100) 

DIMENSION  A (100)  ,B ( 1 01 ) ,C ( 100 ) . D ( 100 ) . E ( 100 ) , F ( 100 ) 

DIMENSION  XU(  100)  .YU(100)  tU(100)  , P(200) 

COMMON/MAI NWJ/PF (3 ,25) ,NPF , STEPS 
COMMON  /MA  INSL/STFT(7)  ♦ LAMSL 

COMMON/HBDB  MB /XI  (100)  , BETA ( 1 00 ) , UE ( 100 ) ,FW ( 100 ) , CPH I ( 100  ), RWC ( 1 00 ) 
S,CF(  100)  ,XC(  100)  »YC ( 100)  ,X(100)  ,CHORD,EPS,TVC,REC,NT,  AK,  Hi,  ETAET, 
$LAM,NPRNT,JMAX, JMI N, JJET, JSEPP 
COMMON/HBDB  BO /XX  (400)  ,YY(400  ) 

EQUIVALENCE  ( U( 1 ) ,UE ( 1 ) ) , ( XU ( 1 ) ,X C ( 1 ) ) , ( YU ( 1 ) , YC( 1 ) ) 

CALCULATE  POTENTIAL  FLOW 

READ  908, CL, ALPHA, RECtNEPS 

READ  909, AK>H i'i»  NT»NB  PTS  .KUTTA .NAEQN 

READ  911.  XDJE  T,S LOTT  ,UTSIG,UTSI  N,  STE  PS , NW  J IT  , LAMS  L 

CHORD  = 14.625 

THICK  = .2 

CAMS  = .05 

C INPUT  BODY  COORDINATES 

CALL  BODY  ( CHORD  ,THI  C K ,C  A MB  , NBPTS  , BT  , RLE,  RT  E ) 

C INPUT  SHOULD  BE  COMPLETE 

- NPRNT  = 0 

NLE  = NBPTS/2  ♦ 1 
PI  = 3.141593 
TOPI  = 6.283185 
RADEG  = 360. /TOPI 
B 2 = 2 • *B  T 

C CALCULATE  PSI  AND  OMEGA 

DO  1 J=l, NBPTS 

HXY  = 1.  - ( XX( J) /B2) **2  - (YY( J)/B2)**2 
SR  = ( HXY**2  ♦ ( YY(J)/BT)**2)**.5 

IF  (HXY  .LT.  0.  .AND.  SR  .LT.  ABS(HXY))  SR  = ABS(HXY) 

SO  = ( .5*(HXY+SR) ) **.5 
J)  * ARSIN(SO) 

0.  .AND.  J .LE.  NLE)  ZO(J)  = PI  - ZO(J) 

0.  .AND.  J .GT.  NLE)  ZO(J)  = PI  ♦ ZO(J) 

0.  .AND.  J .GT.  NLE)  ZO(J)  = TOPI  - ZO(J) 


Z0( 

IF 

IF 

IF 


(XX(J)  .LT. 
(XX(J)  .LT. 
(XX(J)  .GE. 


SHS  = . 5*(  SR-H  XY) 


IF 

SH 

IF 

IF 

ZP( 

IF 


SHS  .LT. 

« SHS**.  5 
(J  .LT.  NLE 
(J  .GT.  NLE 
J)  = A LQG(  SH  ♦ 
(NPRNT  .EQ.  0) 


0.)  SHS  = 0. 


AND.  YY ( J)  .LT.  0.)  SH  = -SH 
AND.  YY ( J)  .GT.  0.)  SH  * -SH 
(l.  + SHS)  **.5) 

GO  TO  3 
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DO  2 J=1,NBPTS 

2 PRINT  901 ,J,XX( J) ,YY( J) ,ZPI J) ,ZOI J) 

3 DPHI  = TOPI /NEPS 
DO  4 J=1 » NE  P S 
P(J)  = (J-.5)*DPHI 

4 PI NEPS+J)  = I J-l J *DPHI 
CALCULATE  EPSILON 
DO  5 J = 1 » NE  P S 
LL  = 1 

IF  (J  .GT.  1)  LL  = LlJ-1) 

DO  5 JJ=LLt NBPTS 
IF  (Z0(  JJ)  .GT.  P(  J)  ) GO  TO  5 
L(  J)  = JJ 

5 CONTINUE 
DO  6 J=1 »NE  P S 

6 F ( J ) = 0. 

JL  = LINEPS) 

JM  = 1 
JU  = L C 1 1 

D 1 = TOPI  - Z0(  JL) 

D2  = ZO(JU) 

D3  = D1  + 02 

Dll)  = -Z  PC  JL)  *D2/D1  /D3  «■ 

DO  7 J=JL. NBPTS 

7 FI  1)  = F( 1)  * IP  I J) 

DO  8 J=1,JU 

8 FI  1)  = FI  1)  ♦ 


- o 


3 


ft  ‘ / 


l1' 


r/ 


ZPIJM)  *U./D1-1./D2)  + ZPIJU  )*D1/D3/D2 


ZPI  J) 

FIl)  = F I 1 ) /( NBPTS- JL+ JU+l ) 


DO  10  J=2  fNEPS 


JU 

= 

LI  J) 

JL 

S 

LI J-l) 

JM 

= 

IJU  ♦ 

JL)  /2 

D 1 

s 

ZOI  JM) 

- ZO I JL) 

02 

X 

ZOI  JU) 

- ZO  I JM) 

D 3 

s 

ZOI  JU) 

- ZOI  JL) 

?P'J 


f 


lU) 


DIJ)  * -ZPI JL) *D2/D1/D3  ♦ ZP I JM)  * 11 . / D1  -1./ D2 ) ♦ IP  I JU  )*D1/D3/D2 
DO  9 JJ=JL.JU 
9 FI  J)  * Ft  J)  ♦ ZPI  JJ) 

10  FI  J)  = FI  J)  /I  JU-JL+l) 

SO  = PI/I2*NEPS) 

JM  * NE  PS/2  - l 
DO  11  JS1 t JM 

SR  = SINI  $0*(2*J+l)l/SINIS0*(2*J-i)  ) 

11  Ct  J)  = ALOG  I ABSI  SR)  ) 

DO  13  J=1 fNEPS 
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» JU  - NEPS 


yt>  ’ 


E(  J)  = DPHI  *D(  J) 

00  12  J J* 1 » JM 
JU  * J ♦ JJ 
IF  (JU  .GT.  NEPS)  JU 
JL  = NEPS  + J - JJ 
IF  (JL  .GT.  NEPS)  JL  = JL  - NEPS 

12  E(J)  = E(J)  ♦ (F(JU)-F(  JL))*C(JJ) 

13  E( J)  = E( J) /PI 
IF  ( NPRNT  .EQ.  0)  GO  TO  16 
DO  1 A J=1 »NEPS 

1A  PRINT  902 ,J,F( J)  ,D( J)  ,E< J) 

00  15  J = 1 ,JM 

15  PRINT  903  «J  »C  ( J) 

CALCULATE  PHI 

16  00  17  J=1 » NE  PS 
P(J)  = P(  NEP  S+ J)  - E(J) 

A(  J)  = OPHI  *<  J-l ) + E(  J) 

IF  { A ( J ) .LT.  0.)  A(  J)  = A(J)+TOPI 

IF  ( A ( J ) .GT.  TOPI)  A ( J)  = A ( J ) -T  OPI 

IF  (P(J)  .LT.  0.)  P(J)  = TOPI  ♦ P ( J ) 

IF  (P(J)  .GT.  TOPI)  P(J)  = P(J)  - TOPI 

17  CONTINUE 

CALCULATE  EPSILON  PRIME 

F ( 1 ) = • 1 *(  - 2.  *E  ( NE  PS-1  )-E(NEPS)+E(2)+2.*E(3)  )/DPHI 

F ( 2 ) = .1*(-2.*E(NEPS)-E(1)*E(3)+2.*E(A)  J/DPHI 

F(NEPS-l)  = .!*(  -2.  *E  (NEPS-3)-E(NEPS-2M-E(NEPS  )*2.*E(  1)  )/OPHI 

F(NEPS)  = • 1 *(  -2  . *E  ( NE  PS-2  ) -E  t NEPS-l  ) + E ( 1 ) +2  .*  E ( 2 ) ) / DP  HI 

LL  = NEPS-2 

DO  18  J=3,LL 

18  F(J)  = .l*(-2.  *E(J-2)-E(J-l  )■*•£(  J+l)-*-2.*E(J+2))/DPHI 
00  20  J=1  »NEPS 

LL  = 1 

DO  19  JJ=LL»NBPTS 

IF  (ZO(JJ)  .GT.  A(  J) ) GO  TO  20 

19  L(  J)  = JJ 

20  CCJ)  = 0. 

IF  (L<  1)  .LT.  L(  NEPS)  ) LL*i 

IF  (L(l)  .LT.  L(NEPS))  GO  TO  22 

DO  21  J=2,NEPS 

IF  (L(  J)  .LT.  L(  J-l)  ) LL=J 

IF  ( L( J)  .LT.  L ( J— 1 ) ) GO  TO  22 

21  CONTINUE 

22  JM  = LL  - 1 

IF  (JM  .EQ.  0)  JM  = NEPS 
JL  = L( JM) 
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23 


24 


25 


26 


27 

28 


29 

30 


JU  = NBPTS 

DO  23  J=JL,JU 

C(  LL)  = C(  LL)  «■  ZP(  J) 

JU  = L(  LL  ) 

DO  24  J=1,JU 
C ( LL ) = C(LL)  + ZP(J) 

C(LL)  = C(LL)/(NBPTS-JL+JU+l) 
AO  = C ( LL ) *DPHI 
26  J=1 »NEPS 


DO 

IF 

JU 

JM 

IF 

JL 

DO 


LL)  GO  TO  26 


0)  JH  = NEPS 


(J  .EQ. 

* L ( J ) 

= J-I 
(JM  .EQ. 

= L ( JM) 

25  J J = JL  » JU 
C(J)  = C(J)  ♦ ZP(JJ) 

C(  J)  = Cl  J)  /(  JU-JL+1) 

AO  = AO  ♦ C ( J)  *DPHI 

CONTINUE 

AO  = AO/TOPI 

IF  (NPRNT  .EQ.  0)  GO  TO  28 

PRINT  903  * LL«AO 

DO  27  J=1 *NE PS 

SO  = PINEPS+J)  *RADEG 

SR  = P(  J)  *RADEG 

PRINT  901,  J,  SO,SR,C(  J)  ,F(  J) 

ELK  = EXPIAO) 

FIND  FRONT  STAGNATION  POINT 

ALF  = ALPHA 

ALPHA  = ALPHA/RADEG 

IF  (KUTTA  .EQ.  1)  CL  = 8. /CHORD*P  I *BT*ELK*S  I N<  E ( 1 ) +ALPHA  1 
AR  = CL/8.*CH0RD/BT/ELK/PI 
SR  = ARSIN(-AR) 

PH  I ST  = ALPHA  ♦ PI  - SR 

PRINT  900 

PRINT  907 

L ( 1 ) * 1 

DO  30  J=2 .NEPS 

LL  = L(J-l) 

DO  29  JJ=LL, NBPTS 

IF  (ZO(JJ)  .GT.  PINEPS+JI)  GOTO  30 
L( J)  = JJ 
CONTINUE 
LL  = 1 

IF  ( PHI  ST  .GT.  PI)  LL  = NEPS/2-1 
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GO  TO  32 

( <E( JSAVE+i)-E( JSAVE)  )/DPHI  )* 


GO  TO  3* 


DO  31  J=L  L » NEP  S 
SO=P(NEPS+J)-E( J) 

IF  ( SO  .GT.  PHI  ST) 

31  JSAVE  = J 

32  ESTAG  = E(JSAVE)  ♦ 

$(  PHI ST-P( JSAVE) ) 

WSTAG  = PHI  ST  + ESTAG 
LL  = 1 

IF  (WSTAG  .GT.  PI)  LL=NLE 
DO  33  J=LL  t NBPTS 
IF  (ZO(J)  .GT.  WSTAG) 

33  JM  = J 

34  PSIST  = ZP(JM)  ♦ ( ( ZP(  JM+1) -ZP<  JM)  )/  (ZO(  JM+1  )-ZO(JM  ) ) )* 

S(  WSTAG-ZO(  JM) ) 

EP  = EXPtPSI  ST) 

EM  = EXP(-PSIST) 

XSTAG  = B T*CQS(  WSTAG  )*(  EP*E  M) 

YSTAG  = BT*SIN(  WSTAG) -MEP-EM) 

CALCULATE  U AND  CP 

PRINT  904,CL,ALF  ?XSTAGtYSTAG 

PRINT  905 

DO  36  J=1 *NE PS 

IF  ( J .NE . 51)  GO  TO  35 

PRINT  900 

PRINT  907 

PRINT  904,CLtALF , XSTAG, YSTAG 
PRINT  905 

35  SP  = P(  J) 

LL  = L ( J) 

SO  = DPHI  *(  J-l) 

PSI  = ZP(  LL)  ♦ ( CZP< LL+l) -ZPCLL)  ) / ( Z0( LL+1 )-ZOCLL ) ) )*<SD-ZOCLL)  ) 
EP  = E XP(  PSI  ) 

EM  = E XP( -P  SI  ) 

TS  = SI  N(  SO) 

TC  = COS(  SO) 

SHS  = • 5*(  EP-E  M) 

SR  * ((1.  + D(  J)  **2)  *(  SHS**2  ♦ TS**2  ) ) **.  5 

IF  (J  .EQ.  1 .AND.  KUTTA  .EQ.  1)  SR  = l.E-8 

SO  = SIN(  SP-ALPHA)  ♦ AR 

U(J)  = ELK*ABS< S0*(1.-FC J) ) )/SR 

PC  J)  = l.  - U(  J)  **2 

SP  = SP  *RADEG 

XU(J)  = 8 T*TC*(  E P+EM) 

YU(  J ) = BT*TS*(  EP-EM) 

36  PRINT  906  » J » SP  f XU(  J)  »YU(  J)  »UC  J)  » P(  J) 
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c 

c 

c 

c 

c 


BYE,  TEDDY 

CHANGE  TO  BOUNDARY  LAYER  COORDINATES 


37 


38 


39 


AO 


Ch 


IF  ( NAEQN  .EQ.  0)  GO  TO  A3 
JM  = NE  P S/2 
B ( l)  = .3A785A85 
B ( 2 ) = . 6 521A515 
8 ( 3)  = B ( 2) 

B(  A)  = B ( I ) 

C( 1)  = -. 86113631 
C(  2)  = 339981  OA  " 

C( 3)  = -C(2I 

ci  a)  = -cm 

DO  A 2 J=1 , JM 
A(J)  = 0. 

A ( JM+  J ) = 0. 

IF  (J  .GT.  1 .AND.  J . LT.  JM)  GO  TO  39 
XLL  = XU(  J) 

XtJL  = XU(J+1) 

00  37  J J=1 , A 

CALL  SLOPE  ( XLL,XUL,C(  JJ)  ,YPU,YPL) 

SR  = • 5/YPU 

A(  J)  = A(  J)  + B ( JJ)  *SR 

A(  J ) = ABS(  .5*(  XUL-XLL)  *A(  J)  ♦ YU(J+1)  - YU(J)) 

XLL  = XU(  JM+  J) 

JL  = JM4-J+1 

IF  (JL  .GT.  NEPS)  JL=1 
XUL  = XU(JL) 

00  38  J J=1 » A „ . 

CALL  SLOPE  ( XLL, XUL, C(  JJ)  ,YPU,YPL) 

SR  = .5  / YP  L 

A(JM+J)  = AIJM+J)  + B ( J J)  *SR 

A ( JM*- J ) = AB  S(  • 5 *(  XUL-XLL)  *A  ( JM+  J)  «-  YU(JL)  - YUCJM+J)) 
GO  TO  A2 
XLL  » XU(J) 

XUL  ■ XU(  J+  1 ) 

00  AO  J J=1 » A 

CALL  SLOPE  (XLL, XUL, C(  JJ)  ,YPU,YPL) 

SR  = ( 1.  + YPU**2)  *♦.  5 

A(  J)  = A(  J)  ♦ B(  J J)  *SR 

A ( J } * ABS(  . 5*(  XUL-XLL)  *A  ( J)  ) 

XLL  * XU(  JM+  J) 

XUL  = XU(JM+jm 


lS 


(Ua 
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41 

42 

43 


44 

45 


46 

47 


DO  41  J J = 1 » 4 

CALL  SLOPE  I XLL»XUL,CI  JJ)  ,YPU,YPL) 

SR  = I i.  + YPL**2)  **.5 

AIJM+J)  = A ( JM+  J ) + BIJJ)*SR 

A ( JM+  J ) = ABSI.5*!  XUL-XLL)  *A<  JM+J)) 

CONTINUE 

GO  TO  52 

LL  = NEP  S/2 

L( l)  = 1 

DO  45  J =2  » L L 

SO  = XU(  J f 

JL  = L(J-l) 

DO  44  JJ=JL,NLE 

IF  ( XX(  JJ)  .LT.  SO)  GO  TO  45 

L( J)  = JJ 

CONTINUE 

LM  = NEP  S/2 

LL  = LM  + 1 

L(  LL)  = NLE 

DO  47  J=LL,NEPS 

SO  = XUIJ) 

JL  = L(J-l) 

DO  46  JJ=JL,NBPTS 
IF  (XX(JJ)  .GT.  SO)  GO  TO  47 
LIJ)  = JJ 
CONTINUE 
LL  = NEPS/2-2 
= NEPS/24-3 
= NEPS-2 
51  J=1,NEPS 

(J  .GT.  3 .AND.  J .LT.  LL)  GO  TO  49 
LM  .ANO.  J .LT.  LU)  GO  TO  49 


48 


LM 

LU 

DO 

IF 

IF 

SO 

JL 

JM 

IF 

JU 

IF 

DO 

SO 

SR 

SO 

JM 

IF 


NEPS)  JM  = 1 


(J  .GT. 

= 0. 

= LIJ) 

= J+l 
(J  .EQ. 

= L(JM) 

(J  .EQ. 

48  J J= JL 

= SO  ( XXI  JJ*1)-XX<  JJ)  ) / CYY  ( JJ*l)-YY(Jjn 
= SO/IJU-JL+l) 

* I l.«-SR**2)  **.5 
= J+l 

I JM  .GT.  NEPS)  JM=1 


NEPS) 
, JU 


JU  = NBPTS-1 


A(  J ) = AB  SI  S0*(  YUI  JM)  -YUC  J)  ) ) 
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GO  TO  51 

49  SO  = 0. 

JL  * L { J ) 

JU  = L(  J*  1 ) 

00  50  JJ=JL,JU 

50  SO  = SO  + ( YY( JJ+1)-YY( J J)  ) / (XX  ( J JU  ) -XX  ( J J ) ) 

SR  = SO/(JU-JLM) 

SO  = ( 1.  + SR**2)**.5 

A(J)  = AB  S(  S0*(  XU(  J*1)-XU(  J)  ) ) 

51  CONTINUE 

52  IF  ( WSTAG  .GE.  PI)  GO  TO  54 
JL  = 2 

JU  = NEP  S/2 
DO  53  J=JLf  JU 
JSAVE  = J 

IF  (XUU+l)  .LT.  XSTAG  .AND.  XU(J)  .GE.  XSTAG)  GO  TO  56 

53  CONTINUE 

54  JL  = NEPS/2+1 
JU  = NEPS-1 
DO  55  J=JL,JU 
JSAVE  = J 

IF  (XU(J«-1)  .GT.  XSTAG  .AND.  XU(J)  .LE.  XSTAG)  GO  TO  56 

55  CONTINUE 

C REINITIAL  U AND  P 

56  P SR  * ALPHA  ♦ ARSIN(-AR) 

IF  { P SR  .LT.  0.)  PSR  = TOPl  + PSR 

DO  57  J=1,NEPS 

0(  J)  = DPHI *( J-l ) - E(  J) 

IF  <B(J)  .LT.  0.)  B(J)  = TOPI  + BIJ) 

57  CONTINUE 

IF  (Bll)  .LT.  B(NEPS))  JM  = NEPS 
IF  (B(l)  .LT.  B(NEPS))  GO  TO  59 
JU  « NEPS-1 
DO  58  J=1,JU 

IF  ( B ( J+l  ) .LT.  B(  J)  ) JM  = J 

IF  ( B ( J+l ) .LT.  B ( J)  ) GO  TO  59 

58  CONTINUE 

59  JU  = JM+1 

IF  (JU  .GT.  NEPS)  JU  = 1 
SR  = B(  JU)  ♦ TOPI 
SO  = BUM) 

IF  (PSR  .GE.  SO  .AND.  PSR  . LT.  SR)  JREAR  = JM 

IF  (PSR  .GE.  SO  .AND.  PSR  . LT.  SR)  GO  TO  61 

DO  60  J=1 t NEP S 
IF  (J  .EQ.  JM)  GO  TO  60 


139 


60 

61 


62 


63 


64 


JU=1 

.AND. 

.AND. 


JU  = J+l 

IF  (J  .EQ.  NEPS) 

IF  (PSR  .GE.  B(J) 

IF  (PSR  .GE.  B ( J ) 
CONTINUE 

IF  (JREAR  .LT.  JSAVEI 
IF  (JREAR  .GE.  JSAVE) 
JBOT  = NE  PS- JT OP 
00  62  J = 1 » J TOP 
JM  = J SA  VE  + 1 - J 
IF  ( JM  .LE.  0) 

B(  J)  = U(  JM) 

C(  J)  = P(  JM) 

D(  J j = XU(  JM) 

FIJI  ~ YU(  JM) 

DO  63  J=1 » JBOT 
JL  = J TOP*  J 
JM  = JSAVE  + J 
IF  (JM  .GT.  NEPS) 


PSR 

PSR 


LT. 

LT. 


B( JU) ) 
B ( JU ) ) 


JREAR  = J 
GO  TO  61 


JT  OP  = 
JT  OP  = 


JM  = NEPS+JM 


JSAVE-  JREAR 
JSAVE  + NEPS -JREAR 


B(  JL)  = U( JM) 
C(JL)  = P ( JM) 

D(  JL  ) = XU(  JM) 

F(  JL)  = YU(  JM) 

01  = XU(  JSAVE)  - 
XU(  JSAVE  + 1) 
01+02 

XSTAG  - XU( 
XSTAG  - XU( 


XU( 

02 

03 

04 

05 

06 
SR 
SR 
JU 

IF  (JSAVE 
X(  l ) = SR 
JU  = JTOP-l 
DO  64  J=i  .JU 
JM  = JSAVE-J 
IF  (JM  .LE.  0) 

X(  J+ll  = XI  J) 
XtJTOP  + l)  = A( 

JU  = JBOT-l 
DO  65  J=1,JU 
JL  = JTOP  + J 
JM  = JSAVE  +J 
IF  (JM  .GT.  NEPS) 


JM  = JM-NEPS 


XU( JSAVE-1 ) 
- XU(  JSAVE) 


+ ( A(  JSAVE) +A(  JSAVE-1)  ) *04/03*0  5/02 


JSAVE-i) 

JSA VE) 

XSTAG  - XU(  JSA  VE+ 1 ) 

-A(  JSAVE-i)  *04/01*06/02 
SR  - A(  JSAVE-1) 

NEPS/2+1 

.EQ.  JU)  SR  * A ( JU  ) * (XST AG-XU  ( JU  ) ) / ( XU  ( JU  + i )-XU(  J U ) ) 


JM  = NEPS+JM 
+ A(  JM) 

JSAVE) -SR 


JM  = JM  - NEPS 
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AND  BETA 
1 ) *U(1)  /CHORD 

5 *X(  JTOP+1 ) *U(  JTOP«-l)  / CHORD 


65  X(  JL  + 1 ) = X(  JL)  * A ( JM) 

DO  66  J = 1 » NE  PS 
U(  J)  = B { J) 

P(  J)  = Cl  J) 

XU(  J ) = D ( J ) 

66  YU( J ) = F ( J) 

CALCULATE  XI 
XI  (1)  = . 5*X( 

XI  ( J TOP  + 1 ) = 

JU  = J TOP-1 
DO  67  J=2  » JU 

XI  = X<  J-l)  /CHORD 

X2  = X(  J)  /CHORD 

X3  = X(  J+l)  /CHORD 

Cl  = U(  J-l) /<  X2-X1) /(X3-X1) 

C 2 = -U( J J /< X2-X1) /(X3-X2) 

C 3 = U(  J + l) /(  X3-X1)/(X3-X2) 

Dl  = C1+C2+C3 

02  = -Cl*(  X2+X3)  - C2  *(  XI +X3  ) - C3*(X1+X2) 

D3  = Ci*X2*X3  *•  C2*X1*X3  + C3*X1*X2 

67  XI  ( J ) = XI  ( J-l ) «-  (X2-X1)*<D1*X1*X2*.5*D2*(X1+X2)*D3) 

$ *D1*<  X2-X1)  **3/3. 

XI(JTOP)  = XI ( JTOP-1 ) + tX3-X2)*(Dl*X2*X3*.5*D2*(X2+X3)+D3) 
$ ♦ D1  *(  X3-X2 ) **3/3. 

JL  = JTOP+2 

JU  = NE  PS-1 

DO  68  J=JL.JU 

XI  = X( J-l) /CHORD 

X2  = X(  J)  /CHORD 

X3  = X(  J+l)  /CHORD 

Cl  = U(  J-l)  /(  X2-X1 ) /(X3-XI ) 

C2  = -U(  J)/(  X2-X1) /(X3-X2) 

C 3 = U(  J + 1)/(X3-X1) /(X3-X2) 

D 1 = C1+C2+C3 

D2  = -Cl*(  X2+X3)  - C2  *(  XI + X3  ) - C3*(X1+X2) 

D3  = C 1*X2*X3  ♦ C2*X1*X3  * C3*X1*X2 

68  XI  f J ) = XKJ-ll  ♦ <X2-X1)*(DI*X1*X2  + .5*D2*(X1«-X2)+D3) 

$ ♦ 01*(  X2-X1)  **3/3. 

XHNEPS)  = XI(NEPS-l)  + < X3-X2 ) * ( D1*X2*X 3+ .5*02* < X2 *X3 )+D3 ) 
$ * Dl*(  X3-X2)  **3/3. 

DO  69  J=1 tNEPS 

69  B(  J ) = X(J) /CHORD 

Dl  = ( 2.*B(  1)*U<1)+2.*B(2)*UI2)-B(1)*U<2)-B(2)*U(1) )/ 

$IB< 1 )**2+B( 2)**2-B(l) *B(2) ) 

BETA(l)  = XI(l)*Ol/U(l)  **2 


i 
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D 2 = 3.  *(  U(  1 ) *B  ( l)  + U(2)*B(2)+U{3)*B(3))-(U{l)+U(2)+U(3l)*(B(  1) 
$ + B ( 2 ) + B(  3 ) ) 

03  = 3.*(B(  1)**2  + B(2)**2  + B(3  ) **2 ) -< B (1) + B< 2 ) +B(3> )** 2 
01  = 2.  *02/D3 

BETA ( 2)  = XI (2) *D1/U<2) **2 

01  = ( 2.*B( JTOP+i) *U( JTOP+i) +2.*8(JTOP+2)*U ( JTOP+2 )-B ( JTOP +1 )* 
$U(  JTOP+2)  -8  ( JTOP+2)  *U(  JTOP+1  >)/( B ( JT OP+i ) **2 +B( JTOP+2 )** 2- 
$B  ( J TOP  + 1 )*B( JTOP+2)) 

BE  TA ( J TOP  + 1 ) = XI ( JTOP+1) *01 /U  ( JT  0P+ 1 ) **  2 
JL  = JTOP+1 
JM  = JTOP+2 
JU  = JTOP+3 

D2  = 3.  *(  U(  JL)  *B(JL)  + U(JM)  *B(JM)+U(JU)*B(JU)  ) 

$ -C  U(  JU  + UI  JM)+U(  JU)  ) *(B(  JL)  + B(JM)+B(JU)  ) 

03  = 3.*(B(  JL)**2+B(  JM)**2+BtJU)**2)-(B(JL)  + B(JM)+B(JU))**2 
01  = 2. *02/03 

BETA!  JTOP  + 2)  = XI  ( JTOP+2  ) *01 /U  ( JTOP+2)**2 
JU  = NEPS-2 
JL  = JTOP-1 
JM  = JTOP+3 
DO  73  J=3  » JU 
IF  ( J .GE.  JL 
DO  70  J J = 1 » 9 
70  C<JJ)  = 0. 

LL  = J-2 
LU  = J+2 

HXY  = .25*<B(  LU)  -BILL)  ) 

XI  = BILL) 

DO  71  JJJ  = LL  ? LU 
DX  » ( B(  JJJI-Xl ) /HXY 
PI  = 1 5*D  X 

P2  = l.-2.*DX+.5*0X**2 
= C(  1)  + P2**2 
= C(  2) 

* C ( 3 ) 

* C(  4) 

= C(  5) 

= C C 61 
= Cl  7) 

= C(  8) 

C ( 91 


AND.  J .LT.  JM)  GO  TO  72 


C<  1) 
C<  2) 
C<  3) 
C(  A) 
C(  5) 
C<  6) 
C(  7) 
C(  8) 
71  C(9) 


+ 

♦ 

+ 

+ 

♦ 

+ 

♦ 

+ 


Pl**2 

1. 

PI  *P2 

P2 

Pi 

U(  JJJ)  *P2 
U(  JJJ) *P1 
U(  JJJ) 


D 1 = C(  7)  - C(  5) *C ( 9)  /C  ( 3 ) 
D2  = C ( 2 ) - C(6)**2/C(3) 

D3  * C(4J  - C(5)*C(6)/C(3) 
D4  = C(  8)  - C(6)*C(9)/C(3) 
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non 


80 


81 

82 


83 


D5  = C(l)  - C<  5)  **2  /C  ( 3 ) 

06  = D 5*0  2 - 03**2 
Cl  = ( D 1*D2-D3  *D4)  /06 
C2  = (04*D5-D1*D3)/D6 

DX  = (C1*CCB(  J)-Xl)/HXY-2.)-.5*C2)/HXY 
2 « *XI  ( J)  *D  X / U I J ) **2 
, JL)  GO  TO  73 

B<  J) -XI)  /HXY-2.  ) -•  5 *C2  ) / HXY 
♦XI  ( J)  *DX/U(  J)  **2 
B(  J+l  )-Xl) /HXY-2.  )-.5*C2)  /HXY 
2.*XI(  J+l)  *DX/U(  J+l ) **2 

8(NEPS-l)-Xl)/HXY-2.)  -,5*C2)/HXY 
) = 2. *XI ( NEPS-1) *DX/U(NEPS-1 )**2 
B( NEP S) -XI ) /HXY-2. )-.5*C2) /HXY 
2.  *XI  ( NEPS)  *DX/U  ( NEPS)  **2 


BE  TA<  J I = 

72  IF  (J  .NE. 
DX  = C C 1 *(  C 
BETAC  J)  = 2 
DX  = CC1*(  ( 
BE  TA  ( J + l ) = 

73  CONTINUE 
DX  = (Ci*(( 
BETA(  NEPS-1 
DX  = (Cl*(t 
BETA(NEPS) 
DO  74  J=1 f 5 
JJ  = JTOP+J 
IF  (BETA(J) 
IF  (BETACJJ 

74  CONTINUE 

J = JTOP+1 
JJ  = JTOP+2 
IF  (BETA!  1 ) 
IF  ( BETA( J) 


.GT.  1.)  BETA  ( J)  = 1. 
) .GT.  1.)  BETA(JJ)  = 


LT.  BET  A( 2 ) ) BETA(l)  = l.+X 1 1 1 * C BET A(2)-l.)/X( 2) 

. LT.  BETA! JJ) ) BETA! J)  = 1 .+X ( J )* C BET A( J J )-l . )/X ( JJ » 


CALCULATE  BOUNDARY  LAYER  FLOW 

DO  80  J = 1 ,NEPS 
FW(  J ) = 0. 

CPHI(J)  = 0. 

RWC(J)  = 1. 

DO  81  J=1 » JTOP 

IF  (XCCJ)  .GT.  XCJET)  GO  TO  82 
JJET  = J 

NPF  * JTOP- J JE  T+ 1 
DO  83  J=1 »NPF 
JJ  = JJET-l+J 
PF ( 1,  J)  = X(  JJ)  /CHORD 
PF  ( 2 1 J)  = U(JJ) 

PF(  3»  J)  = BETA!  JJ)*U(JJ)  **2/ <2.  *XI<  J J ) ) 

XJET  = X(  JJET)  + ( X(  JJET+1)-X(  JJET)  )*(  XCJET -XC  (JJ  ET  ))  / 
$ C XCl  JJET+1)-XC<  JJET)  ) 

XJET  * XJET  /CHORD 
EPS  = 0. 
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c 

c 

c 


90 

900 

901 

902 

903 

904 

905 

906 

907 

908 

909 

911 

912 


TVC  = 0. 

ETAET  = 100. 
NPRNT  ■ 0 
LAM  = 1 
JMAX  = NEPS 
JMIN  = JTOP+1 
CALL  SLAYER 
LAM  = 1 
JMAX  = J TOP 
JMIN  = 1 
CALL  SLAYER 


•O' 


-v ' 


CALCULATE  WALL  JE T F LOW 


STFT(l)  = 
STFT(  2)  = 
STFT(  3 ) = 
STFT(4)  = 
STFT(  5 ) = 
STFT( 6)  = 
STFT(  7)  = 

DO  90  J=1 
STF  T(  6 ) = 
PRINT  912, 
CALL  WALJET 
CONTINUE 
FORMAT  I 1H1 ) 
FORMAT  (I 
FORMAT 
FORMAT 
FORMAT 


SLOTT 

REC 

PF( 2,1  ) 
i.-PF(  2,1)  **2 
XI  ( JJET) 
UTSIG-UTSI  N 
XJET 
,N WJI  T 

STFT(6H-UTSI  N 
J 


6 »4E20. 8) 

{ I 6,3E20.  8) 

( I6,E20.8) 

( 6 X,4HC  L =,F6.3,6X  ,7H  ALPHA  =,F6.2,8H  DEGREES  , 6X,  7HXSTAG  *, 
$E16.8,6X,7HYSTAG  =,E16.8) 

FORMAT  ( / /7H  PO I NT  ,7X  ,3HPHI  ,14X,1HX  , 19X  , 

( I 6 , 5X  ,F  7 .2  ,4  E20 .8  ) 

( 1 OOX ,14HP0TE NTI  AL  FLOW) 

(F6.3,F6.2,E10.3,I3) 

(F9.6,F6.3,I4,I4,I2 ,1  2» 

{ F7.  3,F8.5,F6.3,F7.4,F6.2 ,13,  I2» 

(1H1.19H  WALL  JET  I TERATI  ON,  14) 


FORMAT 

FORMAT 

FORMAT 

FORMA  T 

FORMAT 

FORMAT 

STOP 

END 


1 H Y , 1 9 X , 1HU, 19X, 2HCP  ) 


SUBROUTINE  BODY  (CH ORD  tTH  I CK ,C AMB»  N8 PTS»  BT » RL  E,  RT  E ) 
COMMON/HBDBBO/XX(400)  ,YY(40D) 

KIND  AIRFOIL 
Cl  = .5+CHORD 
C 2 = C 1 *THI C K 
N1  = NBPTS/2 


DTHET  = 3.141593/N1 
A = .065 

RTE  = C1*THICK**2 
JTE  = 0 
XXI  1)  = CUA 
XX( N 1+ 1 I = -Cl  + A 
YY(  1)  = 0. 

YY( N1 *1 ) = 0. 

DO  1 J=2  » N1 
JJ  = NBPTS+2-J 
G = ( J- 1 ) *D  THE  T 
XX(J)  = A+C  1*C0S(G) 

YY(  J)  = C2*$IN(G) 

XX(JJ)  = XX(J) 

YY(  JJ)  = - YY(  J) 

IF  (JTE  .EQ.  0 .AND.  XX(J)  .LT.  6.815)  JTE=J-i 
CONTINUE 

IF  ( JTE  .LT.  2)  GO  TO  3 
RTE  = .5525 
CHORD  = CHORD-. 01 
XX( 1)  = 7.3675 
DTHET  = 1.570796/1  JTE-1 ) 

DO  2 J=2 • JTE 
JJ  « NBPTS+2- J 
G = <J-1)*DTHET 
XX(J)  = 6.  81 5+R  TE  *C  0S(  G ) 

YY ( J ) = RTE*SIN(G) 

XX(JJ)  = XX(J) 

YY(JJ)  = - YY(  J) 

RLE  = C1*THICK**2/12. 

RTE  = RTE/12. 

CHORD  = CHORD/12. 

DO  4 J=1 »NBPTS 
XX(J)  = XX(J)/12. 

YY{  J ) - YY(  J)  /12. 

BT  = .25*(CH0R0-. 5*C RLE+RTE)  ) 

RETURN 
END 


TS  3V 


cr  c.  < ^ 


kxcjo-  „ 7a)'L{£, 


1 4-*  (f7$ 


C Uott-v.  - . o l 


(4. 


t 
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SUBROUTINE  SLOPE  ( XL L »XUL  *X  A ,Y PU  , YPL ) 
X = .5*1  XUL-XLL)  *XA4-.5*(XUL«-XLL) 

C KIND  AIRFOIL 

X - i 2 . *X 

IF  (X  IGE.  6.815)  GO  TO  I 
XS  = X-.065 

SR  = ABSC  2.  13 89 062 5-. 04 *XS* *2 ) 

SR  = SR**.  5 
YPL  = .04*XS/SR 
YPU  = -YPL 
RE  TURN 

1 XS  * X-6.815 

SR  = ABSI  .3 052 5625-X S**2  ) 

SR  = SR**. 5 
YPL  = XS/SR 
YPU  = -YPL 
RETURN 
END 


vw-^ 


l W 


V 


Jr 


o 


SUBROUTINE  BLAYER 

DIMENSION  E I 5,241)  .B(2A1  ,A)  ,Y(2A1)  ,C  (81)  ,G(2,  241) 

COMMON/BL  SLBL/EDVI5  ,2A1 ) , NBL.NEDGE 

COMMON /HBDBMB /XI ( 1 00 ) ,BE TA ( 1 00 ) ,U E ( l 00) , FW ( 100 ),  CPHI ( 100 ) ,RrfC( 100) 
$ »CF(  100)  » XC  ( 100)  »YC ( 100)  »X(100I  ,CHORD,EPS,TVC,REC,NT, AK, HI,  ETAET, 


SLAM,NPRNT,JMAX, JMI N,JJET, JSEPP 
COMMON /HBDBEB /A ( 5 ,241 ) ,0(60)  .ETAI2A1 ) 

,A)  ,F(2,2A1) ,H(2A1) 


( 5,2A1 

COMMON/HBDBMP/U(2A1 
COMMON  PI  4,241) 

DOUBLE  PRECISION  A ,E  ,P  ,U  , B,Y  , ET  A , H,  C , D,  F ,G 

PRINT  900 

JCOUNT  = 0 

DU)  = EPS 

D(  2)  = TVC 

Dll)  = REC 

D(  A)  = 0(3) **.5 

D<  5)  = Oi 3) **.25 

D(  A 01  = 0. 

DIA7)  = CHORO 
CALL  PROFLE 
JU  = 1 
JL  * 1 
D(  8)  = 0. 

0(9)  = 0. 

D(  10)  = 0. 

D(  11)  = 0. 

D(  16)  = 0. 

D(  17)  = 1. 

D(  39)  = 0. 

GO  TO  2 r/'J 

JU  = JMA  X 
JL  = JMIN 
DO  99  JPT  = JL , JU 
JCOUNT  = JCOUNT*  1 
IF  (JCOUNT  .NE.  1) 

IF  (LAM  ,EQ.  0)  GO 
N = 2A1 
LU  = N-2 
D(  6)  = 6. 

D(  7)  = D(  6)  /(  N-l ) 

DO  3 J=l , N 
H ( J » = 0(7) 

E TA(  J)  = ( J-ll  *D(7) 

C(  1)  = 1. 

GO  TO  6 


GO 

TO 


TO 

A 


10 


1A7 


4 N = NT 
LU  = N-2 
Cl  1)  = AK 
D(  6)  = ETAET 
D( 7)  = HI 
E TA  (1)  = 0, 

ETAI2)  = D(  7 ) 

HI  1)  = HI 
H ( 2 1 = Cl  1) *Hi 
C(  2)  = C(  1 ) 

00  5 J-3.N 
C ( 2 ) = C I 2 I *C  1 1 ) 

ETA(J)  = DI7)*IC(2)-1.)/(C(1)-1.) 

Cl  1)  *HI J-l) 

Cl  1) **2 
C<  11  *C(  2) 

C(  l)  *C(  3) 

C(1)*C(4) 

C(  1 1 *C(  51 
l.+CI 1) 

CI7H-C12) 

Cl  8)  +CI  3) 

= 1./ICI7)  *CI8)*C(9)  I 
= -1.  /(C(7)*C(8)*C(3)) 

= 1.  /<C(  7 1 **2*C  I 5)  I 
= -l./(C(7)*C(8)*CI6)) 

= 1 ./I  C I 7)  *C  I 8)  *CI9)  *C  (6  ) ) 

= C « 1 0 1 /C  C 14) 

= C ( 1 1 1 /C ( 14) 

= C ( 1 2 ) /C  ( 14) 

= C(13)/C(14) 

= Cl  18)*CI9)  /C(8) 

* :(17)+C(9)/C(7I 
= C I 1 61  *C  I 9) 

= C(3)*CI7)*C<8)  + 

$ C I 7)  *C  I 81  *CI  9) 

C(  23)  = Cl  18)*CI7)  *C(8)*C(9)/C(22) 

: I 17) *C(1) *C(8) *C(9) /C(22) 
C(16)*C(2)*CI7)*C(9)  /C  (22 ) 

C(  15)*C(3) *C(7)*C(8) /C(22) 
6.*IC<7)+CI1)-C(2I ) 
2.*C(l)*(C<7)-C<7) *C(1)-CI2) ) 

-Cl  3 ) *C I 7) 

6.*<C(7)+C<i)-CI2)*C<7) ) 

2.*CI  1)*CI7)*(1.-C(1)*CI7)-C(2) ) 


5 HI  J) 

6 C(  2) 

C(  3) 

C ( 4) 

C(  5) 

C(  6) 

C<  7) 

C(  8)  = 
C(  9)  > 
C(  10) 
C<  11) 
C(  12) 
Cl  13) 
C(  14) 
C(  15) 
C<  16) 
C(  17) 
C(  18) 
C<  19) 
C(  20) 
C(  21) 
Cl  22) 


Cl  24) 
C<  25) 
C I 26 ) 
Cl  27) 
Cl  28) 
Cl  29) 
Cl  30) 
Cl  31) 


C 12 ) *C (7 )*C  (9  ) + 


Cl  1)*C(8)*C(9)  ♦ 
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t 


C C 32)  = -C(3)*C(7>**2 

C(  33)  = 6.*(C(7)+C(l)-C(2)-C(2)*C(7)  ) 

C( 34 ) = 2.*C( 1) * 

$ ( C ( 3 ) *C  ( 7)  - C ( 2 ) *C  ( 7 ) - CC2)  - C(i)*C<7)**2  - C<1)*C(7)  ♦ C(7)) 
C( 35)  = Cl 3) *C< 7) *(C (2)  ♦ Clt)*C(7)  - C<7)  - 1.) 

C ( 36 j = 6.*(C(7)  - C ( 2 ) - C(2)*C(7J) 

C ( 37)  = 2.*C(2)*C(7)  *(C(2)-C<7)-1.) 

C(  38)  = C<4)*C<7)**2 

C(  39)  = 6«*C(1)*(1.-C(1)-C(1)*C(7)) 

Cl  40)  = 2.*C13)*IC11)*C17)-C(7)-1.) 

C( 41)  = Cl  5)*Cl  7) 

C ( 42 ) = 6.*IC(  8)+C(l) *C(7)+C(2) ) 

C C 43)  = 2.*C(1)*(C(7)*C18)+C(8)*C11)+C17)*C<2)  ) 

C<  44)  = C13)*C17)*C18) 

C ( 45 ) = 6.  *1  C(8)  +C  1 1 ) *C(7)+C(2)-C(3)  ) 

Cl 46)  = 2.  *C  ( 1 ) * 

$ i C ( 7)  *Cl  8 ) + C(8)*C(l)  - C ( 8 ) *C  ( 2)  + C<7)*C12)  - C ( 7 )*C  ( 3)-Cl  4)  ) 

C ( 47)  = -C13)*IC13)  *C17)+C12)*C18)*C11)*C<7)*C18)-C17)*C18)) 

C ( 48)  = 6.*IC18)+CU)*C17)-C<3) ) 

C ( 49)  = 2.*C11)*IC17)*CIB)-C12)*C(8)-C17)*C(3)) 

C(  50)  = -C14)*C17)  *C18) 

C(  51)  = 6. *1  Cl  8)  ♦ C (2  ) -C  ( 3 ) ) 

C ( 52)  = 2.*C12)*IC(8)-C18)*C11)-C13)  ) 

C( 53)  = -C( 5)*Cl 8) 

C ( 54)  = 6.*Cll)  *IC17)+C11)-C12) ) 

C(  55)  = 2.*C(3)*lC(7)-C(i)*C(7)-C(2)  ) 

C(  56)  = -CC6)  *Cl  7) 

C ( 57)  = C(10)*IC17)*C(8)*C19)  ♦ C(8)*C(9)  ♦ C17)*C19)  ♦ C17)*C18)) 
C(  58)  = C(  11)*C17)  *C18)  *C(9) 

C112)*C18)*C(9) 

Cl  60)  = C ( 13) *C ( 7) *C  1 9) 

Cl  14)  *C  l 7)  *C18) 

-C( 10) *C (3) ♦C  f 7 1 *C18) 

Clll)  *C(2)  *l-C  1 ) *C(7 )*C (8 ) ♦ C (7)*C( 8 ) ♦ C ( 8 ) ♦ C€  7 ) ) 

C(  12)  *C(  2)  *C  (7)  *C  18 ) 

C113)*C(2)*C(8) 

CI14)*C(2)  *CC 7) 

C ( 3) *C ( 7) 

Cl 10)*C( 2) 

Clii)*Cli)*Cl7) 

C< 12)*(C (2)+C<lt *Ct7)-C(7)-l.) 

-Cl 1 3 ) *C  l 7 ) 

C( 72)  = -C( 10) *C ( 3) *C(  7) 

C(  73)  = -Cl  1 1 ) *C  ( 2)  *C  l 8) 

C ( 74)  = -Cl  12) *C<1) *C17)*C<8) 


C(  59) 
C(  60) 
C(  61) 
Cl  62) 
C(  63) 
C(  64) 
C(  65) 
C(  66) 
C(  67) 
Cl  68) 
C<  69) 
C(  70) 
C(  71  j 
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C C 75  I = -C(  13 ) *(  C ( 3) *C(7)+C(2)*C(8)+C(i)*C(7)*C(8)-C(7)*C(8) ) 

C( 76)  = C< 14) *C( 7) ♦CCS) 

C C 77)  = C(10)*C(3)*C(7)  *C(8) 

C ( 78)  = C( 11) *C( 2) *C(7) *C(9) 

C ( 79)  = C(12)*C(1)*C(8)*C(9) 

C(  80)  = C(  13)*C(  7)  *C(8)  *CC9) 

C ( 81)  = C(  14)*(C(3)*C(7)  *C(B)  + C (2  ) *C  (7  ) *C  ( 9 ) ♦ C ( 1 )*C ( 8 )*C ( 9) 
$ C ( 7 ) +C  ( 8 ) *C  ( 9 ) ) 

DO  7 J = 1 y N 
G(  l, J)  = 0. 

7 G(  2.  J)  = 0. 

IF  (LAM  .EQ.  0)  GO  TO  9 
DO  8 J=1 t N 
F(  1,  J)  = 0. 

8 F ( 2 1 J ) = 0. 

GO  TO  10 

9 CONTINUE 

CALL  PTIRB  <N) 

10  IF  ( JU  -EQ.  1)  SO  TO  16 
LU  = JPT 
LM  = JPT-1 
LL  = JPT-2 

IF  (JPT  ,GT.  JL  .AND.  JCOUNT  .EQ.  1)  GO  TO  14 
IF  (LL  .GE.  JL)  GO  TO  13 
IF  LM  -EQ.  JL)  GO  TO  12 
D(  8)  = I./XKJL) 

D(  9)  = -D ( 8) 

D(  10)  = 0. 

GO  TO  15 

12  D(  8 ) = 1 • /(  XI  ( JL+ll-XI  ( JL)  ) + 1-  / XI  ( JL+1 ) 

D(  9)  = -XII  JL  + 1 ) /XI (JL) /(XI (JL+1)“XI(JL) ) 

0(10)  = ( XI  ( JL+1  )-XI  ( JL)  ) /XI  ( JL) /XI  ( JL-H  ) 

GO  TO  15 

13  D(  8)  = l./(  XI(LU)-XI  (LM))+i./(XI(LU)-XI(LL)) 

0(9)  = -(  XI  (LU)-XI  <LL))/(XI(  LU)-XI(LM))  / (XI  (LM)-XI(LL)  ) 

DC  10)  = ( XI(LU)-XICLM)  )/(XICLM)-XI  ILL))  / (XI(LU)-XI(LL)) 

GO  TO  15 

14  DC  S)  = 1.  /(XI  ( LU)-XI  ( LM)  ) 

D(  9)  = -DC  8) 

D(  10)  » 0. 

15  D( 11)  = XI ( JPT) 

D(  12)  = ( 2- *D( 11 ) ) **.5 
D( 13)  = UE( JPT) 

D( 14)  = SPHI ( JPT) 

D(  15)  = RWC(  JPT) 


0(16)  = 2.*DC 1)*D(  2) *D(12)*D(14) /D(l 5 ) **2/ D( 13 )/ 0(4 ) 

0(  17)  = BE  TA ( JPT ) 

0( 39)  = D( 15) 

IF  (D(  1)  .LT.  .5)  D(  39)  = 1. 

0(39)  = -0  ( 12) /D(4)/D(13)/D(39) 

LU  = N-2 

16  00  40  JI T = 1 ,20 
TOL  * .0001 

IF  (JIT  .GE.  10)  TOL  = .001 

IF  (LAM  .EQ.  0 .AND.  JCOUNT  .EQ.  1)  TOL  = .015 
JPR  * JIT 
00  17  J=1.N 

E( 1, J)  = 1.  + D( 16) *ETA( J) 

E ( 2 , J ) = D(  16)  ♦ P(1,J)  ♦ ETA(J) 

E ( 3 , J ) = -D( 17) *CP(2 ,J)+2.)  - 2.*Ddl)*D(8)*(P(2,J)+l.J 
E ( 4,  J ) = 2.  *D(  11  )*P(3»J)*D(8) 

17  E(  5,  J)  = -<E(3,J)*P(2»  J)+E(4.  J)*P(1,J))  ♦ 2 .*D(  1 1 )*  ( ( P ( 2,  J ) ♦ 1 . )* 
*(  D(  9)*F(  2 , J)+D<  10)  *G  (2  , J)  ) -P(3  , J)*(0(9)*F(  1,  J)+D(  10)*G(  1,  J ) ) ) 

IF  (LAM  .EQ.  1)  GO  TO  22  ^ 

DO  18  J*1,N 
A ( 1,  JI  = P(  1 , J) 

A ( 4,  J ) = P ( 2 , J)  * ' 

A ( 2 , J ) = P(  3 , J) 

18  A(  3, J ) = P( 4,J) 

CALL  EODIE  (N) 

00  19  J=1,N 
ED V(  5,  J)  = 1 • + A ( 1 , J) 

E(  1 , J 1 = E(  1 » J)  *(  l.  + AC  1 • J)  ) 

19  E(  2 , J ) = EI2,J)*D(16)*A(1  , J)  ♦ A(  2 , J ) * (1  . + D (16  )*ET  A(  J ) ) 

22  00  23  J=1,N 

E(  5,  J)  = c(  5 » J ) - ( E (1  ,J)  *P  (4  » J)  «•  E (2  ,J)*P(3»J  ) ) 

E(  2 , J ) * E(  2,J)/E(1  , J) 

E( 3, J)  = E(  3, J) /E(l ,J) 


T 


0 / L,  ■ 1 


.4- 


23 


A( 1,2)  = 
A( 2 ,2 ) = 
A(  3,2)  = 
A(  1 ,N  i = 
A ( 2»N ) = 
A(3,N)  = 
A( 4,N) 


E(4,J)/E(1,J) 
E ( 5 , J ) /E<1  ,J) 


C(  19) 

C(  20) 

C ( 2 1 ) 

C(  23) 

CC24I 
C(  25) 

C ( 26) 

0(  18)  = E(2,N-1)  *H(N-4) 

D(  19)  = E(3 ,N-i) *H(N-4) **2 

0(20)  = C(42)+D( 18)*C(43)+D(19)*C(44) 


. 
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A ( 1 » N- 1 ) = C ( 1 8 ) *(C(45)*-D(18)*C(46)+D(19  )*C (47 ) ) / D(  20  ) 
S «-E(4,N-l)*H(N-4)  **3/C(  14)/D(20  ) 

A(  2 »N-  1)  = C<17)*(C(48)«-D(18)*C(49)*D(19)*C(50))/D(20) 
A ( 3 » N- 1 ) = C(16)  *(C(51)+D(18)*C(52)*D(19)*C(53)I/D(20) 
A(4,N-1)  = C<15)*(C(54)+D(18)*C(55)*0(19)*C(56)I/D(20) 
A(  5 »N-1 ) * E <5 ,N-1) *H(N-4)**3/C(14)/D(20 ) 

DO  24  J=3,LU 

D(  18)  = E(2,  J)  *HC  J-2) 

D(  19)  = E (3,  J)  *H  ( J-2) **2 

D(  20)  = C( 27)+D(  18)*C(28)«-D(19I*C(29) 

D(  21 ) = H( J-2) **3 

A ( 1 , J)  = C<  18)*(C(30H-D(18)*C(31)+D(19)*C(32))/D(20) 

J)  * C(  17)  *(C(33H-D(18)*C(34)+DC19)*C(35))/D(20) 
♦E(  4,  J)*D(  21)/C(14)  /D(20) 

J)  = C(16)*(C(36)+D(18)*C(37)+0(19)*C(38))/D(20) 

J)  = C(15)*(C(39)+D(1B)*C(40)+D(19)*C(41))/D(20) 

J)  = E( 5, J) *0(21 ) /C(14) /D(20) 

1)  = A < 3 * 2 ) 

2 = A(  2 *2 ) 

3)  » A( 1 ,2) 

4)  = 1. 

1)  = A(  3,3)/U(2,l) 

1)  * A(  2 *3  ) -B  ( 3 *1 ) *U(2  *2  ) 

2)  = A(l»3) -8(3*1) *U( 2 *3 ) 

3)  = l.-B(3,l) 

2)  = A(4,4) /U(2,l) 

1)  = (A(3»4)-B(4f2)*lM2f2))/U(3tl) 

1)  = A(2 ,4)-B(4,2) *U(2,3)-B(4,1)*U(3,2) 

2)  = A(  1,4)-B(4,2)-B(4  ,1  )*U(3,3) 

3)  = 1. 


A ( 2 

$ 

A ( 3 
A ( 4 

24  A(  5 
U(  2 
U(  2 
U<  2 
U(  2 
B(  3 
U<  3 
U(  3 
U<  3 
B(  4 
B(  4 
U(  4 
U(  4 
U<  4 
DO 

B ( J 
B(  J 
U(  J 
U(  J 

25  U(J 


f J I — 1 • 

25  J=5 • LU 

,2)  = A(4,J) /U(  J-2  ,1 ) 

,1)  = (A(3,J)-B(  J,2)  *U<  J-2  ,2)  )/U(  J-l,  1) 

,1)  = A(2*J)-8(J,2)*U(  J-2  ,3)-B(J*l)*U(J-l»2) 

,2)  = A( 1 * J) -8 ( J * 1 ) 

,3)  = 1. 

= A(  4»N-1 ) /U(N-4  ,11 

= (A(3,N-1)-B(  N-l  ,3  ) *U  ( N-4 ,2))/U(N-3,  1) 

* ( A(2,N-l)-B(N-l  ,3  ) -B  (N-l  ,2 ) *U  ( N-3 ,2)  )/U  (N-2,  1) 
= A(  1 ,N-1 1 -B (N-l ,2)-B(N-l,l)*U(N-2,2) 

= 1 • — B ( N-l ,1) 


B(  N-l, 3) 

B ( N- 1 ,2 ) 

BC N-l  ,1) 

U(  N- 1 ♦ 1 j 
U(  N-  1 *2) 
B(N,4)  = A(  4 
B(  N *3) 

B(  N , 2) 

B(  N , 1 ) 


A(4,N)/U(N-4,1) 

(A(3,N)-B(N»4) *U( N-4  »2 ) )/U(N-3*l) 

( A(2,N)-B(N,4)-B(N,3)*U(N-3,2) )/U(N-2, 1) 
(A( 1,N)-B( N*3) -8(N»2) *U ( N-2 ,2) )/U(N-l» 1) 
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U(N,1)  = l.-B(N,2)-B(N,l)*U(N-i,2) 

Y(  2)  = 0. 

Y(  3)  = A( 5,3) 

0D  26  J=4,LU 

26  Y(J)  = 4(5, J)  - B(  J,2)  *Y<  J-2  ) - B( J, 1 )*Y (J-l ) 

Y(N-l)  = A(  5.N-1  )-B(N-l  .3)  *Y  (N-4) -B(N-1,2)*Y  (N-3  ) -B ( N-  1.  1)*  Y ( N-  2 ) 
Y(N)  = -B(N,4)*Y(N-4)  -B(N,3)  *Y(N-3)-B(N,  2)*Y  (N-2)-B(N,  1)*Y(N-1) 

A ( 1 ,N  ) = Y(  N)  /U(  N »l ) 

A ( 1 , N- 1 ) = ( Y(  N-l ) -U  ( N-i  ,2)  *A(  l ,N)  )/U(N-l,  1 ) 

00  27  J=4,LU 

1 = LU+4-J 

27  A(  1 , 1 ) = (Y(I)-A(l»I+2) -U (I  ,2) *A (1,1+1)  )/U(I,l) 

A(l, 3)  = ( Y(3)-U(3,3)  *A(1  ,5) -U(3,2)*A(1,4)  )/U(3,l) 

A(  1,2)  = ( Y(2)-A(l,5)-U(2,3)  *A(1  ,4 ) -U  ( 2 , 2 ) * A ( 1,  3 > >/U  ( 2,  1) 

A(  1,1)  = 0. 

DO  28  J=2  ,3 
K = J-l 

A(J,1)  = -(C(57)*A(K,1)+C(58)*A(K,2)+C(59)*A(K»3)+C(60)*A(K,4) 

S *C(  61)*A{  K,5)  ) /H(l) 

A(  J ,2  ) = (C(62)*A(K,1)+C(63)*A(K,2)+C(64)*A(K»3)*-C(65)*A(K,4) 

$ +C(  66) *A(  K,5) ) /H(l) 

A(J.N-l)  = C(3) *(C(72) *A(K,N-4)+C(73)*A(K,N-3)+C(74)*A(K,N-2) 

$ +C( 75)*A(K,N-1 ) ♦C( 76 ) *A ( K, N) )/H(N-4) 

A ( J , N ) = C ( 3 ) *( C ( 77)  *A(K,  N-4  )«-C(78)*A(K,N-3)+C(79)*A(K,N-2) 

$ *C( 80)*A(K,N-i )+C(81) *A( K, N) )/H(N-4) 

00  28  J J = 3 » LU 

28  A(  J , J J ) = C ( 67)  *(C(681  *A(  K,  JJ-2  )+C(69)*A(K.  JJ-1 1«-C(70)*A(K,  JJ  > 

$ +0(71)  *A(  K,JJd)-C(14)  * A ( K , J J<-2 ) ) / H ( J J -2  ) 

DO  29  J = 1,N 

29  4(  4,  J)  = E(5,J)-E(2,  J)  *A(3  , J ) -E  ( 3 , J ) *A(  2 , J ) -E  (4,  J)*A(  l,  J ) 

IF  (LAM  .EQ.  1)  GO  TO  31 

IF  (JIT  .GT.  1)  GO  TO  31 
00  30  J=1 ,N 
A ( 1 , J ) = . 5*A(  1 , J) 

A(  2,  J)  = • 5 *A  ( 2 » J) 

A(  3 , J)  = • 5*A ( 3 , J) 

30  A ( 4 , J ) = .5*A(4,  J) 

31  IF  (JU  .EQ.  1)  GO  TO  32 
0 ( 22)  = F W(  JPT) 

P(  1,1)  = D(22) 

32  00  33  J=1,N 

P(  1,J)  * P(  l,J)+A(l  ,J) 

P(  2 , J ) * P(  2 * J)  + A(2  , J) 

P(  3 , J ) = P(  3 » J)+A(  3 , J) 

33  P(  4,J)  = P(  4 , J)  + A(4  , J) 


I 
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IF 

( P ( 3 » 1 ) .LE. 

0.  ) 

JSE  PP= 

: JPT 

IF 

( P( 3.1)  .LE. 

0.  ) 

PRI  NT 

906 , JPT 

IF 

( P ( 3,1  ) .LE. 

0.  ) 

GO  TO 

104 

DO 

34  J=2 . N 

IF 

( P(  2 , J ) .GT. 

0.  ) 

GO  TO 

35 

IF 

( P( 2, J-l ) .GT 

. P ( 2 , J)  ) 

GO  TO  35 

34  JSAVE  = J 
GO  TO  37 

35  JLOW  = J SAVE* l 
DO  36  J= JLOW  ,N 

P(  1 , J)  = P(  1 .JSAVE) 

P ( 2 » J ) = 0. 

P( 3.J)  = 0. 

36  P<  4 »J ) = 0. 

37  CONTINUE 

0(23)  = (2.*A(3,1))/(2.*P(3,1)-A(3,l)> 

DFW  = 0(23) 

OF  W = ABSIDFW) 

JPRNT  = JIT 
DOS  = 0. 

IF  (LAM  .6Q.  I)  GO  TO  21 
IF  (0(2)  .GT.  .5)  GO  TO  38 
0(24)  = D(39)*(P(l  ,N)-P(1  ,1)  ) 

GO  TO  20 

38  D( 24)  = 0. 

00  39  J=2  »N 

39  0(  24)  = 0(24)  * . 5*D(  39)  *H  ( J-l ) * ( P(2  , J)  ♦ P (2  , J- 1 ) )/ 
$( l.+.5*D(  16) *(ETA(  J)+ETA( J-I ) ) ) 

20  D( 40)  = (0(24) -D (40) )/0(24) 

DOS  = 0(40) 

OOS  = ABS(ODS) 

0(40)  = 0(24) 

21  IF  (DFW  .LT.  TOL  . ANO.  DOS  . LT.  TOL)  GO  TO  41 

40  CONTINUE 

41  IF  (JU  .EQ.  1)  GO  TO  48 
D(  24)  =1. 

IF  (0(1)  .GT.  .5)  D ( 24)  = RWC(JPT) 

0(24)  = 2. *0(24)  *P(3  ,1 ) /D(4) /D(  12  ) 

CF( JPT)  * 0(24) 

DO  42  J=1,N 
JM  = J 

IF  ( P(  2 . J)  .GT.  -.005)  GO  TO  43 

42  CONTINUE 

43  IF  (0(1)  .GT.  .5)  GO  TO  44 
0(  41)  = -D(  39)  *E  TA(  JM) 


45  0(42) 


46  0 ( 43) 


D(  42)  = D(39)*(P(1 »JM)-P(1 »1)) 

GO  TO  46 
44  D(  41)  = 0. 

0(42)  = 0. 

DO  45  J=2,JM 

0(43)  = (ETA(J)-ETA(J-1))/(1.  + .5*D(16)*(ETA(J)+ETA(J-1)))**.5 
D(  41)  = 0 ( 41 ) +0 ( 43) 

= 0(42)  + .5*(P!2,J)*P(2*J-1))*D(43) 

D( 41)  = -D( 39) *0(41) 

0(  42)  = D ( 39) *0(42) 

= 0. 

DO  47  J=2»JM 

47  0(  43)  = 0(43)  ♦ . 5 *(  ET  A ( J)  -ETA  ( J-l ) ) * ( P ( 2,  J M-P  ( 2»  J- 1 ) )*  ( 1 . + . 5* 

% (P( 2,  J)  + P( 2, J-l) ) ) /(I .*.5*D(16)*(ETA( J )+ET A(J-l) ) )**.5 

D(  43 ) = 0 ( 3 9)  *D ( 43  ) 

• 0( 44)  = D( 42 ) /D( 43) 

D(  45)  = 0(3)  *0(13)  *0(42) 

0(46)  * 0(3) *0(13) *0(43) 

IF  ( NPRNT  .EQ.  1)  PRINT  900 

PRINT  901  » JPT,X(  JPT)  ,XI  (JPT)  , BET  A ( J PT ) ,U  E ( JPT  ) » X C(  JPT  ),  YC  ( JP  T ) 
PI  = 12. *D(  411*0(47) 

P 2 = 1 2 . *0 ( 42)  *0(47) 

P 3 = 12. *0(43) *0(47) 

P4  = CF(JPT) 

P 5 = D( 44) 

PRINT  902  »P1 »P2  »P3  *P4  *P5 
PI  = 0(3) 

P2  = 0(45) 

P3  = 0(46) 

PRINT  903  .JPRtDFW»DDS»Pl  »P2  t P3 
PI  = P( 3,1) 

IF  (NPRNT  .EQ.  0)  PRINT  909,JM,P1 
GO  TO  50 

48  PRINT  904  » JPR  ,0F  M 
DO  49  J=1  t N 

F ( 1, J)  = P( 1 ,J) 

49  F ( 2 * J ) = P(  2,J) 

50  IF  (NPRNT  .EQ.  0)  GO  TO  52 
PRINT  905 
00  51  J=1 »N  » 1 0 

ETA(J) 


PI 
P 2 
P3 
P4 
P 5 


P(  If J) 
P(2,  J) 
P( 3, J) 
P(4,  J) 


E TA  ( J) 
1. 


l/l 

Ln 


JET 


51 

52 


53 

54 

55 


PRINT  908,P1  ,P2,P3,P4,P5 
IF  ( JU  .EQ.  1)  SO  TO  99 
IF  (JPT  .EQ.  JJET)  GO  TO  101 
DO  53  J=i,N 
G(  It  J)  = F(  1 , J) 

G(2,J)  * F(  2 t J) 

F ( 1,J)  = PC  1 , J) 

F ( 2,J)  = PI  2 t J) 

IF  (LAM  .EQ.  1)  GO  TO  55 
GO  TO  99 

IF  (PC  3,1 ) .GT.  .21  GO  TO  60 
IF  I D(  46)  .LT.  125.)  GO  TO  99 
IF  (JPT  .EQ.  JU)  GO  TO  99 
D( 48)  = -.001371*0(46)+. 0817 

0(49)  = D(47)*(UE(  JPT+1)-UE(  JPT ) ) / (X  ( JPT  + 1 ) -X  ( JPT  ) > 

D ( 49 ) = D(3)*0(43)**2*0(49) 

IF  (0(49)  .LT.  D ( 48)  ) GO  TO  60 
LAM  = 0 
JCOUNT  = 0 
PRINT  911 
GO  TO  99 
D( 48)  = 640. 

IF  (BETA(JPT)  .LT.  0.)  D(48)  = 320. 

IF  (0(46)  .LT.  D (48 ) ) GO  TO  99 
LAM  = 0 
JCOUNT  = 0 
PRINT  907, JPT 
CONTINUE 
IF  (JU  .EQ.  1) 

CONTINUE 
GO  TO  104 
00  102  J=1,N 
EOV(  1 ,J) 

ED  V(  2 , J ) 

ED  V(  3 , J ) 

ED  V(  4,  J) 

NBL  = NT 
NEOGE  = JM 
CONTINUE 
FORMAT  (1H1) 

FORMAT  ( //,6H  POI NT, I 3 , IX ,3HX  = , E15 . 8 , 2X , 4HX I = , E14.7. 2X, 6HBETA  =, 
$E12.5,2X,4HUE  = ,E  14.  7 ,2X  ,4HXC  =, E14.7 ,2X,4HY C =,E14.7) 

902  FORMAT  ( 1 0X*7HDELTA  = ,F8.4,3H  l N,2X,7HDELST  =,F8.4,3H  IN,2X, 

$7H  THE  TA  =,F8.4,3H  IN,2X,4HCF  = , E14.7  ,2X  , 3HH  =,E14.7) 

903  FORMAT  ( 1 OX  ,9H  I TE  RAT  I ON, I 3 ,8  X ,5HDFW  = , El  3.6,  2X,  5 HDDS  =,E13.6,2X, 


60 


99 

100 

101 


102 


104 

900 

901 


GO  TO  1 


E TA  ( J ) 

P ( 1 , J)  + E TA  ( J) 
P( 2 , J) +1 . 

, J) 


= P ( 3 
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904 

905 

906 

907 

908 

909 
oil 
912 


S5HREC  =,E13.6,2X,6HREDS  = ,El  2 . 5 ,2X  ,5  HRET  =,E13.6l 


FORMAT  1 1 8H  STAGNATION  POINT:  ,5X,9HITERATI0N.I3,5X»  5H0FW  =,E13.6) 
FORMAT!  / , 1 1 X , 3HE  TA  ,19X,1HF  ,19X,2HF'  f 17X  ,3HF*  • , 17X,  4HF*  • • ) 

FORMAT  !/,35H  BOUNDARY  LAYER  SEPARATION  AT  POINT, 16) 

I/, 2 OH  TRANSITION  AT  P0INT,IS) 

( 5E2  0.7) 

( 10X  , 7HNE0GE  = ,1  4 ,9X  ,6HFW  • =,EL2.5) 

/54H  'SHORT  BUBBLE'  LAMINAR  SEPARATION:  TRANSITION  ASSUMED) 
( 5E15.7) 


FORMAT 

FORMAT 

FORMAT 

FORMAT! 

FORMAT 

RETURN 

END 
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SUBROUTINE  PROF  IE 
COMMON  P( 4.241) 

DOUBLE  PRECISION  P ,E  TAE  ,E  TAD  ,ET  A , ET  , Cl , C2  » C3  ♦ C4,  C5,  C6,  H 
N = 241 
ETAE  = 6. 

ETAD  = 3. 

Cl  = 4. /ETAD 
C2  = 12./ETAD**2 
C 6 * - . 2*E  TAD 
H = ETAE/(N-1) 

DO  2 J = 1,N 
ETA  = ( J-l ) *H 

IF  (ETA  .GT.  ETAD)  GO  TO  1 
ET  = ETA/ETAD 
C 3 = E T **2 
C 4 = ET**3 
C 5 = E T**4 

E TA*( 2*C5+C4-2 .*C3+2. *ET-1. ) 

-C  5+4. *C4-6. *C3+4.*ET-l . 

C1*(-C4+3.*C3-3.*ET+1.  ) 

C2*(-C3+2»  *6 T-l  • ) 

C 6 

0. 

0. 

0. 


P(  1,  J)  = 
P( 2, J)  = 
P( 3, J)  = 
P ( 4 » J ) * 
GO  TO  2 
P(ltJ)  = 
P(2,J)  = 
P( 3, J)  = 
P ( 4 • J ) = 
CONTINUE 
RE  TURN 


» 
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SUBROUTINE  PTURB  ( N) 

DIMENSION  D(  7) 

COMMON  PI  4,241) 

COMMON /HBDBMP/UI 241 ,4)  ,F(2  ,241)  ,H(241) 
DOUBLE  PRECISION  P,U,F,H,D 
JL  = 1 
Dll)  =6. 

D(  2)  = 6.  /240. 

0(3)  = 0. 

DO  3 J=2,N 
0(3) 

JS 


, IN 

D(  3) 


♦ H(  J-l) 


D(  7)  = PC  1 * J-l ) 

IF  ( D(  3)  .GT.  6.  ) GO  TO  4 
DO  1 J J = JL » N 
JM  * JJ 

D(  4)  = ( JJ-1)  *0(2) 

0(5}  = J J *D  ( 2 ) 

IF  (0(3)  .GE.  D ( 4)  .AND.  D(3)  .LT.  D(5))  GO  TO  2 

1 CONTINUE 

2 D(6)  = (D(3)-D(4))/(D(5)-D(4)) 

JL  = JM 

IF  ( JL  .GE.  241)  GO  TO  4 

U(J,1)  = P( l ,JL)+D(6) *(P(1 , JL+1 )-P(l ,JL) ) 

U(  J , 2 ) = P(2,JL)+D(6)*(P(2,JL*-l)-P(2,JL>  ) 

U(  J , 3 ) * P(  3,JL)+D(6)  *(P(3  , JL+1  ) -P(3,JL)  ) 

U(  J, 4)  « P(4,JL)+D(6)*(  P(4  »JL+1)-P(4,JL)  ) 
CONTINUE 
DO  5 J=JS,N 


U(  J,l) 
U(  J,2) 
U(  J,3) 
U(  J,4) 
F ( 1,1) 
F( 2,1) 


D(  7) 

0. 

0. 

0. 

P(  1 ,1) 
P(  2,1  > 


DO  6 J=2,N 

U(  J,l) 
U(  J,2) 


P( 1,J)  = 
P( 2, J)  = 


P(  3 , J } = U(  J ,3 ) 
P(  4 , Jl  = U(  J ,4 j 
F ( 1,J)  = PI  1 , J) 


6 F ( 2» J ) 
RETURN 
END 


P(  2,J) 
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SUBROUTINE  EDO  I E ( N) 

COMMON /HBDBEB/Al 5 ,24 L ) ,D( 60)  , ET A (241 ) 
OOUBUE  PRECISION  A ,0  ,E  TA 
20  OR  AS  WITHOUT  TVC 
0(31)  = 0 { 15) 

IF  (0(1)  .UT.  .5)  0 C 31 ) = 1. 

0(25)  = .16*0(4)  *D(12)/D(31) 

0(26)  = -0(  5)  *(D(  12) /D(  31 ) ) **.5/26. 

0(  27)  = .01 68*0  ( 4)  *D(12)/D(31) 

0(28)  = A(2 ,1 ) 


-.005)  GO  TO  2 


D ( 29)  = A(  1 ,1  ) 

DO  1 J= 1 i N 
D(  30)  = A( 1 , J) 

IF  (A(4,J)  .GT. 

1 JSAVE  = J 

2 JSAVE  = JSAVE  + 1 
00  3 J=1 , N 
0(31)  = A ( 2 , J) 

IF  (0(31)  .UT.  0.)  0(31)  = -0(31) 

3 A(  1,J)  = D{  31) 

0(30)  = (D(29)-D(30)  ) *0(27) 

0(31)  = E TA( JSAVE) 

00  4 J=1 ,N 

0(32)  = D ( 28)  -0  ( 1 7)  *ETA(  J) 

IF  ( D ( 32 ) .UT.  0.)  D( 32  ) = -0(32) 

0(  32)  = D ( 32 ) **. 5 

0(33  = DEXP(0(26) *ETA( J) *0(32) ) 

A(4,J)  = D(  25)  *ETA(  J)  **2*A(1  , J)*(l.-D(33)  )**2 
D(  34)  = 2.*A(2»J)  + ETA(J)*A(3»J) 

IF  ( A( 2«  J ) .UT.  0.)  0(34)  = -D(34) 

A(5,J)  = D(  25)  *ETA(  J)  *11.-0(33)  )*(D(34)*  11.-0(33)  )-2.*ETA(J  )* 
$A(  1 , J ) *0 ( 26)  *0(3 2)  *0(33)  *(  1 . -.  5 *D  ( 17  ) *ET  A(  J )/0  ( 32)**2  ) ) 

0(35)  = 0(30) /( i.+5.5*(ETA(J)/D(31) )**6) 

IF  ( A ( 4 » J ) .GT.  0(35))  GO  TO  5 

4 JSAVE  = J*1 
PRINT  10 

5 0( 36)  = -33. /O (31) 

DO  6 J= JSAVE  ,N 
0(37)  = E TA  ( J ) /O  ( 31 ) 

D(  38)  = 1. *5.5*0(37)  **6 
A ( 4,  Jl  = 0(  30)  /D  (38) 

6 A(5,J)  * 0(36)  *0(37)  **5  *A  ( 4 » J ) / D ( 38  ) 

A(  1,1)  = 0. 

A ( 2,1)  = 0. 

A ( 1,2)  = <A(4,2)+A(4,3))/3. 
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A ( 2,2)  = (A(5,2)+A(5,3) )/3. 

LU  = N-2 
DO  7 J=3 , LU 

A ( 1 , J ) = ,2*(A(4  , J-2H-AI4  , J-i)  + A(4,J)+A(4,  J + 1J+A14, J+2)) 
7 A<  2,  J ) = • 2*(A(5,J-2)-*-A(5,J-l)+A(5»J)  + A(5,J  + l)+A(5»J+2)) 
A ( 1 » N- 1 ) = A ( 1 , N-  2 ) 

A(  1 ,N ) = A ( 1 .N-2) 

A ( 2,N-1)  = A 1 2 ,N-2 ) 

A ( 2,N)  * A ( 2 ,N-2 ) 

10  FORMAT  (19H  EDDIE'S  IN  TROUBLE) 

RE  TURN 
END 
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SUBROUTINE  WALJE  T 

DI  ME  N SI  ON  CPWJ15  0)  »CFWJ(  50)  » UEW  J (50  ) ,XWJ  (50),  DELTA!  50  ) 
DIMENSION  ET1121)  ,E(  4, 4 ,120)  ,Z  (4,121  ),DG(4,12i) 
DIMENSION  A(4,4)  ,B  ( 4 ,4  ) , C (4  ,4  ) , D (4  ,4  ) , Q ( 4 ) 

C GMMON/MAI N WJ/PF (3 ,2  5 ) , NPF  , STE  PS 
COMMON/JE  TMI  X/HSU21 ) ,F(6  ,121  ),S(  18  I 

COMMON /HBDB  JE/ETA(  121)  ,E  PS  (2  ,1 2 1 ) ,G(  4 , 12  I ) , GG (7 , 121  ) 
INPUT  HERE 
CALL  SLOTBL 
JU  « 121 
S(  6)  = S(7) 

S(  8)  = STE P S*S(  l ) 

ET(  1)  = 0. 

ETA(l)  * 0. 

UEP  = S(  4) 

NT  = JU-1 
LU  = JU 
LT  = NT 
MXMN  = 0 
NMA  X = 50 

INPUT  INITIAL  F,G,H,CP  AND  ETA  GRID 


DO  l J=1,JU 
G(  1 , J ) = F(  1 ,J) 
G(  2, J)  = F ( 2 , J) 
G(  3,  J ) = F ( 3 , J) 


4 ; J ) 
F(5,  J) 


DO  2 J=2 , JU 
ETA(  J ) = ETA(  J-l  )+HS(  J) 

ET(J)  = . 5»(ETA( J)+ETA( J-l) ) 

REC  = S(2) 

CS  = S(l) 

CA  = C S**2*REC 
CB  = . 5*C  S*REC 
CPWJ(  1)  = G ( 4 , 1 ) 

UEWJ(l)  = S(  4) 

XSLOT  = S( 6) 

XWJ(  1)  = SC  7) 

DX  = SC  8) 

DXM  = • 5*D  X 

CFWJ(l)  = 2.*S(4)+G(3,1)/(CS*REC) 
DELTA!  1)  = S!  5 ) 

X = XWJ!  1 ) - XSLOT 
CALL  RADI  US  (X,RC) 

GE  = CS/RC 
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DO  3 J=1,JU 

GG( 3» J)  = GE/(1.+GE*ET( J)) 

Cl  = DELTA ( l)  /CS 
JMA  X = 1 
JMIN  = 1 
DO  4 J=2  * NT 

IF  ( G( 2 1 J ) .LT.  G ( 2 » J-l ) .AND.  JMAX  .EQ.  1)  JMAX=J-1 
IF  (G(2,J)  .LE.  G ( 2 » J-l ) .AND.  G f 2 * J I .LE.  G(2,J+1»)  JMIN=J 
IF  (ETA(J)  .GT.  Cl)  GO  TO  5 
JDEL  = J 

IF  (JMAX  .EQ.  1)  JMAX= JDE  L 
IF  (JMAX  .GT.  JDEL)  JMAX  = JDE  L 
DO  6 J=l,4 
B ( 3, J)  = 0. 


NU  = NPF-1 
N 1 = l 

BEGIN  SOLUTION  FOR  POINT  N 

DO  75  N=2  »NMA  X 

PICK  X FOR  STATION  N 

UEPS  = UEP 

JGOB  = 0 

IF  (N  .LE.  3)  JG0B=1 
D XT  = DX 

IF  (N  .EQ.  2)  D XT=C  S/3. 

IF  (N  .EQ.  3)  D XT=C  S 
X = XWJ( N-l ) +D  XT 
DO  7 J=NL  »NU 

IF  (PF( 1 » J)  .LE.  X .AND.  PF(1,J+1)  .GE.  X)  GO  TO  8 
N1  = J 

IF  ( N 1 .GE.  NU)  N1=NU-1 
N2  = Nl+1 
N3  = N2+1 

C A = PF(1,N2)-PF(1,N1) 


C( 2 » J ) 
NL  = 1 


C 2 = -C7*C9/C4/C6 
C 3 = C7*:8/C5/C6 


UEP  * C1*PF(2,N1)+C2*PF(2  t N2  ) *C3  *PF  (2  » N3  ) 

DUEP  = C1*PF(3,N1)  + C2*PF(3,N2)  + C3*PF(3,N3) 

JGOB  = JGOBU 

IF  (JGOB  .EQ.  2)  GO  TO  LO 

DELU  = ABS(  UEP-UEPS) 

IF  (DELU  .LT.  .1)  GO  TO  10 
DXT  = • 1+DX/DE  LU 
IF  (DXT  .LT.  DXM)  OXT=DXM 
X = XWJ( N-l ) + DXT 
GO  TO  9 

10  XWJ(N)  = X 
CO  = 3. 

Cl  = PF(1,N1)+PF(1,N2)+PF(1,N3» 

C2  = PF(  1 ,N1  j**2+PF{  1 .N2)  **2  + PF  (1  , N3 ) **2 
C 3 = PF  ( 3 »N1  J + PF  ( 3 »N2  )*PF  (3  »N3  ) 

C4  = PF(1,N1)*PF(3,NI)  + PF(1.N2)*PF(3,N2)+PF(1,N3)*PF(3,N3) 
ODUEP  = ( C 0*C4-C l *C3 ) / ( C0*C2 -Cl  **2 ) 

X * X-XSLOT 
CALL  RADIUS  (X,RC) 

GE  = CS/RC 

DO  11  J=1 »JU 

GG( 1 ,J)  = i.+GE  *ETA( J) 

GG( 5,J)  * GG(3  * J) 

GG(2,J)  = GE  /(  1 . +GE  *E  TA ( J)  ) 

GG(  3 » J)  = GE/(1.*GE*ET(  J)  ) 

GG  ( 4 y J ) = 2.  *GG(  1 » J)  **2/CA/  (l.+GE*ET  ( J)  ) 

11  CONTINUE 

STORE  N-l  QUANTITIES 
00  12  J=2tJU 

Cl  = . 5*1  G(  1 ,J)+G(  1 , J-l)  ) 

C 2 = . 5*(  G(  2. J)+G(2 , J-l) ) 

C 3 = .5*(G(  3tJ)*G(3  » J-l ) ) 

C4  = .5*(G(4,J)tG(4,J -D) 

F ( 1, J)  = UEHJ(  N-l)  *C1 
F(2,J)  = C2 

F(  3, J)  = C3+GG ( 5 » J) *C2 
F ( 4 »J ) = C4 

12  CONTINUE 

ADJUST  N-l  PROFILE  FOR  INITIAL  GUESS  AT  N 
ETGE  = 1. 2*E  TA  ( JDE  L) 

IF  (ETGE  .GT.  ETA(JU))  ETGE*ETA( JU) 

DO  13  J=JDEl,JU 
LU  = J 

IF  (ETA(J)  .GE.  ETGE)  GO  TO  14 

13  CONTINUE 
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1A  LT  = LU-1 

DELTA!  N)  = DELTA(N-l) 

E 7 = ETGE 

IF  (JMIN  .GT.  I)  E7=l.l  *ETA(  JMI  N) 

Cl  = l.*GE*E7 
USDU  = DUEP/UEP 

C 2 » - • 5*RC**2*DDUE  P*ALOG  (Cl  ) **2 

UWJS  = ( UEP+C2) /Cl 

IF  ( MXMN  .GT.  0)  GO  TO  15 

UEWJ(N)  = UEP 

CPE  = 1 .-  UE  P **2 

GO  TO  18 

15  UE WJ( N)  = UWJS 

Cl  * UE  WJ(  N)  *G(  l » LU) 

DO  16  J= JMA  X,  JU 

IF  I F ( 5 » J ) .GT.  Cl)  GO  TO  17 

16  JM  = J 

17  CPE  = F(6.JM)-UEWJ(N)**2-(RC*DUEP*AL0G(GG(liLU)  )/ GG( l , LU  ) )** 2 

18  Cl  = CPE-G(A.LU) 

DO  19  J=1 *L U 

19  G(4,J)  = G(  4,JH-C1 
JM  = LU+1 

IF  (JM  .GT.  JU)  JM=JU 
DO  20  J=JM,JU 

G(  1 « J ) = G(  1 » J- 1 ) +HS(  J)  *G  (2  » J-l ) 

G( 2,J)  = 1. 

G(  3, J)  = 0. 

20  G( A» J 1 * CPE 
BEGIN  ITERATION 
DO  61  I T=i » 1 0 

L = N 

IF  (N  .EQ.  2 .AND.  IT  .EQ.  1)  L=l 
E A = (G(A»1)-CPWJ(N-1) ) /DXT 

CALL  EDDY  ( C S ,REC  ,UEW  J<  N)  ,EA  ,X , RC , E7  , JMAX  , JM  IN,  J DEL  , LU,  L , IT  ) 
CO  = MXMN 


DO  AO  J=2  » LU 

ci  

C 2 
C 3 : 

CA 
C 5 
C6 
C 7 
C9  i 
CIO 


«-G(l  , J-l)  ) 
♦ G( 2 t J-l ) ) 
+ G(3,  J-l ) ) 
♦G(A, J-l) ) 

*C2 


• 5*(  G(  1,J’ 

. 5*(G( 2, J 
. 5*( G ( 3 » J 

• 5*(  G(  A » J 
C 3+GG(  3,J1 
UE  WJ  ( N ) *C  1-F  ( l , J) 

UE  WJ(  N)  **2  *C2 
HS(  J)  /DXT 

F(  A,J)+(  UEWJ(N-l)  *F  (2»  J)  ) **2-CA-(UEWJ  (N)*C2>**2 


i 

! 

i 
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C12  = 2.*HS(  J)*GG(3,J)*UEWJ(N)**2*C2 
C 1 2 = C12*(1.-.5*CS*PC*ETGE*USDU**2) 

IF  (ETA(J)  .GT.  1 • 2 *E  7 ) C12=0. 

S(  1)  = 5+UE  W J(  N)*C9*(UEWJ(N)*C5HIEWJ(N-i)*F(3,J)) 

S(2)  = C9*(C7-.5*UEWJ(N)  *GG(3,  J ) *C6  ) +UEW J( N) *EPS ( 2,  J)*GG<  2,J  ) 

$ -UEWJ(N)  *EPS<2  , J)  *GG(6.J) 

S( 3)  = -UEWJ(N) ♦t.5*C9*C6+EPS(2 » J) ) 

S(  A)  = • 5*C9 
St  5)  = -C 12 

S<6)  * C9*(C7-.5*UEWJ(N)  *GG(3,  J ) *C6 ) -UEW J ( N ) *EPS ( 2,  J-1)*GG(  2,  J 
$ + UEWJ(N)  *EPS(2.  J-l)*GGt6, J-l) 

St  7)  = UEHJ(N)  *( -.5*C9*C6*-EPS<2  , J-l)  ) 

S(8)  = -G(A,J)+G(A,  J-1)  + C12*C2 

St  9)  = C9*(C10*C6*(UEWJ(N)*C5+UEWJ(N-1)*F(3,  J ) ) ) +U  EWJ  ( N )* 
StEPSt  2,J)*(G(3,J)-GG(2,J)  *G  ( 2 , J)  + GG  (6 , J * * IG  ( 2,  J >-l . I ) 

$-EPS(  2,  J-l)  *(G(3  ,J-1)-GG(2  , J-i  )*G(2,  J-1)+GG(6,  J-1)*(G<  2,  J-l)-  1 
IF  (J  .EQ.  2)  GO  TO  25 

S(  10)  = -G(l,J-l)+G(l,J-2)*-.5*HS(J-l)*(G(2,  J-l  )*G(  2,  J-2) ) 

St  11 ) = -G( 2, J-l )+G(2, J-2)+.5*HS( J-l )*(G(3,J-1 )+G<3,  J-2)) 

25  CONTINUE 

C 1 1 = . 5*H  St  J-l) 

IF  (J  .GT.  2)  GO  TO  31 
DO  26  K=1  »A 


A ( 1 » K ) 
A(  2,K) 
26  A ( A, K ) 
A ( 1,1) 
A ( 2,2) 
A(  A,  A) 
A( 3,1) 
A ( 3,2) 
A( 3,3) 
A ( 3, A) 
A ( A, 2) 


0. 

0. 

0. 

1. 

1. 

-1  . 
S(  1) 
S(  6) 
St  7 1 
St  A) 
S(  5 ) 


« i *»» c i = ai  : 

DO  2 7 K=1  .A 
Ct  3,K)  = At  3,K) 
27  Ct  A,K)  = At  A,K) 
Ct  3,2)  » St  2 ) 

Ct  3,3)  = St  3) 
C(A.A)  = 1. 

Qt  1)  = 0. 

0(2)  = 0. 

Qt  3)  = St  9) 

Qt  A)  = St  8) 

DO  28  K=1,A 


-1) 


.)  ) ) 
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DO  28  KK=  l » 4 

28  D(K,KK)  = A ( K , KK) 

CALL  INVERT  CD) 

DO  29  K=1  ,4 

DO  29  KK=1 » 4 

29  E ( K»KK  , 1 ) = D(  K,1)*C(1,KK)«-DCK,2)*CC2 
$ +DC  K,4)  *C  < 4,KK) 

DO  30  K=l,4 

30  Z ( K » 1 ) = DC K,1 ) *QC I)tO(K,2)*Q(2)+D(K, 
GO  TO  40 

31  DO  32  K=1 ,4 


AC  l.K) 
AC  2»K ) 
AC  4fK) 
B ( 1 , K) 
32  B(  2,K) 
AC  1,1) 
AC  1,2) 
AC  2,2) 
AC  2,3) 
AC  3,1) 
AC  3,2) 
AC  3,3) 
AC  3,4) 
AC  4,2) 
AC  4,4) 
BC  1,1) 
B ( 1,2) 
BC  2,2) 
8 ( 2,3) 


0. 

0. 

0. 

0. 

0. 

1. 

-C  11 

1. 

-Cll 
SC  1) 
SC  6) 
SC  7) 
SC  4) 
SC  5 ) 
-1. 
-1. 
-Cll 
-1. 
-Cll 


D 1 3 I - -LI 

DO  33  K=1 ,4 
CC  3,K ) = AC  3 , K) 
33  CC  4,K)  = AC  4 , K) 
C ( 3,2)  = SC  2 l 


CC  3,3) 
CC  4,4) 


= SC  3) 

= 1. 


H — 1 • 

Q ( 1)  = SC  10) 

OC  2)  * SC  11) 

QC3)  = SC  9) 

QC  4)  = SC  8) 

DO  34  K=1 ,4 
DO  34  KK=1 * 4 

34  DC  K,KK ) = ACK,KK)-BCK,l)*E(l,KK,J-2)- 
$ -BC  K,3)*E(  3,KK,J-2)-BCK,4) *EC4,KK, J 
CALL  INVERT  CD) 


,KK)+DCK,3)*CC3,KK) 

3)*Q(3)+DCK,4)*Q(4) 


BCK , 2 )*  El  2,  KK,  J-2) 
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DO  35  K=l,4 
DO  35  KK=1,4 

35  E ( K » KK » J- 1 ) » D(  K,l)  *C  ( 1 , KK)  «-D  l K ,2 ) *C<2  , KK)  ♦DIK , 3 )*  C ( 3,  KK  ) 

$ +D(  K,4)  *C(4,KK) 

DO  36  K=1 ,4 

36  Q(K>  = Q(K|-BIKfl)*Zll,J-2)-BIKt2)*ZI2,J-2)-BIK,3H'Z(3,J-2) 

$ -B(  K,4)  *Z(  4,J-2) 

DO  37  K=1 ,4 

37  Z(KtJ-I)  « D(K,1)  *Q(l)+D(K,2)*Q<2)+D(K,3)*Q(3)+D{K,4>«‘Q(4) 
40  CONTINUE 

DO  41  K= 1 ,4 


A ( i,K) 
A ( 2, K) 
A ( 3,K) 
A<  4,K) 
3(  i,K) 
B ( 2,K) 
41  Q(  K)  = 
A ( 1,1) 
A( 1,2) 
A ( 2,2) 
A ( 2,3) 
A(  3,2) 
A ( 4,4) 
B ( 1,1) 
B( 1,2) 
B ( 2,2) 


0. 


0. 

0. 

0. 

0. 

0. 

0. 


1. 

-.5*HS(  LU) 

1. 

A ( 1 ,2  ) 

1. 

1. 

-1. 

A ( 1,2) 

B( 2i 3)  = All, 2) 

0(1)  = -G(1,LU)+G(1,LTH-.5*HS(LU)*(G(2,LU)+GI2,LT)) 

Q(  2)  = -G(2,LU)+G(2,LT)+.5*HS(LU)*(G(3,LU)*G(3,LT)) 

DO  42  K=l,4 
DO  42  KK  = 1 ,4 

42  D(K,KK)  = A(K,KK)-B(K,l)*E(l  , KK,  LT  ) -B ( K , 2 ) *E<  2,  KK,  LT  ) 

$ -BI K,3) *E(3,KK,LT) -B( K ,4) *E <4 , KK, LT ) 

CALL  INVERT  ( 0 ) 

DO  43  K = 1 ,4 

43  Q(K)  = Q( K)-B( K,l) *Z(1 ,LT) -B( K, 2 ) *Z<2, LT )-B<K,3)*Z( 3, LT ) 

$ -B(K,4) *Z(4,LT) 

DO  44  K=1 ,4 

44  DG(K.LU)  = D(K,i)*Q(l)  + DIK»2)*Q(2)+D(K,3)*Q(3)+DIK,4)*Q(4) 
DO  45  J J=1 » LT 

J = LU-JJ 
DO  45  K=1  ,4 

45  DG(K,J)  = Z(  K.  J)-E(K,1  , J)  *DG(1  , J+l ) -E< K.  2,  J )*DG<  2,  J +1 ) 

$ -E( K,3, J)*DG(  3 ,J+i)-E  <K,4,  J)*DG<4,  J+l) 


168 


OG I 1,1)  =0. 

DG(  2*1)  = 0. 

DO  48  J=1 ,LU 
DO  A 8 K=l,4 

IF  (IT  .GT.  I)  GO  TO  48 
DG(K,J>  = • 5*DG ( K,  J) 

48  G(K,J)  = G( K, JI+DG( K,J> 

TEST  FOR  SEPARATED  PROFILE 

IF  ( G(  3*1)  .LE.  0.)  PRINT  1 1 4 , N,  I T ,G  (4,  l ) 

IF  <G<  3,1  ) .LE.  0.  ) GO  TO  86 
FIND  JMA  X » JMI  N , JOE  L 
DO  A 9 J=2,JU 
JMA  X = J 

IF  ( G ( 2 , J ) .LT.  G(  2 » J-l ) ) GO  TO  50 
IF  (J  .EQ.  JOEL)  GO  TO  53 
A9  CONTINUE 
GO  TO  53 

50  DO  51  J=JMAX,JU 
JMIN  = J 

IF  ( G ( 2 , J ) .GE.  G(2,J-m  GO  TO  52 

51  CONTINUE 

52  IF  (JMIN  .GE.  JOEL)  GO  TO  55 

53  IF  ( JMA X .GE.  JDEU  JMIN=1 
DO  5A  J=JMIN,LU 

JDEL  = J 

IF  ( G( 2 t J ) .GT.  .99)  GO  TO  57 
5A  CONTINUE 
GO  TO  57 

55  JMIN  = 1 

DO  56  J= JMA  X , JU 
JDEL  = J 

IF  { G(  2 , J ) .LT.  1.011  GO  TO  57 

56  CONTINUE 

57  IF  ( JMA X .GT.  JDEL)  JMAX  = JOE  L 
IF  ( MXMN  .GT.  0)  JMI  N=1 

IF  (JDEL  .EQ.  JU)  GO  TO  60 
IF  (JMIN  .GT.  1)  G (2 , JDEL)  = l. 

G(  1,  JDEL)  = G(  1 ♦ JDEL-1  ) + .5*HS(JDEL)*(G(2,JDEL)*G(2»  JDEL  - 1 ) ) 

ADJUST  SOLUTION  FOR  NEW  DELTA 

JSAVE  = JDEL+l 

DELTA!  N)  = C S*E  TA(  JDEL) 

DO  59  J=JSAVE,JU 

G(  1 , J ) = G(  1,J-1)^HS(  J)*G(2,J-1) 

G(  2, J)  = 1. 

G( 3, J)  = 0. 
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59  CONTINUE 

G(  3»  JOEL ) = (G(2,JDEL+i)-G(2,JDEL-l))/(ETA(JDEL+l)-ETA(JDEL-l)) 

60  G(  3,  JDEL-I)  = (G(2,  JOEL) -G(2  ,JDEL-2)  ) / < ET  A < J DEL  )-ET  A l J DEL-2 ) ) 

IF  (JMAX  .LT.  JOEL)  JMAX  = JMAX-I 

IF  (JMIN  .GT.  I)  JMI  N= JMI  N-l 
C TEST  FOR  CONVERGENCE 

IF  ( I T .EO.  1)  GO  TO  61 
TOL  = .005 

IF  (IT  .GT.  5)  TOL=.  01 

IF  (IT  .GT.  7)  TOL=.  01 5 

DHW  = ABS(DG(3,n/(G(3,l)  + .5*DG(3,im 

OCPW  = ABS(DG(4,1)) 

IF  (DHW  .LT.  TOL  .AND.  DCPW  .LT.  TOL)  GO  TO  62 

61  CONTINUE 

62  CONTINUE 

C OUTPUT  AND  STORAGE 

IF  (N  .LE.  3)  GO  TO  92 
IF  (JMIN  .EQ.  1)  GO  TO  92 

IF  ( G(  2 » JMI  N)  *UE  P .LT.  .98*UWJS)  GO  TO  92 

MXMN  = MXMN+1 

UEWJ(N)  = UEP*G  ( 2 , JMI  N) 

DO  91  J=1 »JU 

F(  5 * J)  = UEP*G(  1 tJ) 

F(  6*  J ) = G(  4 » J)  ( UE  P*G  (2  »J)  ) **2 
IF  (J  .GT.  JMIN)  GO  TO  90 
G(  If  J)  = G(  1 , J)  /G  ( 2 » JMI  N) 

G(  2 » J)  = G(2,J)/G(2,JMIN) 

G(3,J)  = G(  3 » J)  /G(  2 » JMI  N) 

GO  TO  91 

90  G( 1»J)  = G( 1 » J-l )+HS( J) 

G( 2 * J ) = 1. 

G( 3, J)  = 0. 

G(  4,  J)  = G(4,  JMI  N) 

91  CONTINUE 
JOEL  = JMIN 
JMIN  = 1 

92  CONTINUE 
CPWJ(N)  = G ( 4 « 1 ) 

CFWJ(N)  = 2.  *UE  WJ(N)*G(3»l)/(CS*REC) 

UMAX  = G(  2»  JMAX) 

UMIN  = G(  2 » JMI N) 

CPOT  = l.-UEP«*2 

PRINT  100*N  »XWJ(  N)  ,UEWJ(N)  ,CPWJ(N)  , CPE,  CPOT 
PRINT  101  ,CFWJ(  N)  ,X,DELTA(N)  ,UMAX,UMIN 
PRINT  102,1  T,ETA(  JMAX)  , ETA  (JMIN)  , DHW,  DCPW 
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r 


PRINT  103, JOEL, G(3,l) 

NPRNT  = 0 

IF  ( NPRNT  .EQ.  0)  GO  TO  64 
PRINT  104 

C 1 = 2.  *UE  W J(  N)  /(  C S*REC  ) 

DO  63  J = 1 ,LU»3 

C 2 = Cl*EPS(i,J)*(G(3,J)-GG(2,J)*G<2,J)+GG(6,J)*(G(2,  J)-l.)) 

63  PRINT  105,ETA(  J)  ,G(1  ,J)  ,G<2,  J)  ,G(3,Jl  ,G  (4  , J I , EPS  ( 1,  J),C2 

64  CONTINUE 

ADJUST  GRID  FOR  DOWNSTREAM  CALCULATIONS 
IF  (N  .NE.  3)  GO  TO  75 
PRINT  115 
DO  65  J=2  » JU 

IF  (ETA(J)  .GT.  .05)  GO  TO  66 

65  JK  = J + l 

66  CONTINUE 
RS  = 1.055 

GG< 1 » JK-1 ) = E TA ( JK-1 ) 

JEO  = JU 

DO  67  J=JK,JU 

HSU)  = H$(J-1)*RS 

IF  (HSU)  .GT.  1.)  HS(J)=1. 

GG(ltJ)  = GG  ( 1 , J-l ) ♦HSU) 

IF  ( GG(  1 , J ) .GT.  ETA(JDEL)  .AND.  GG(1,J-1)  .LT.  ETA(JDEL))  JEO=J 

67  CONTINUE 
JB  = 1 


DO  70  J- JK, JEO 
DO  68  JJ=JB , NT 


JA  = J J-l 
IF  (ETA(JJ)  .GE.  GG  ( 1 , J)  ) GO  TO  69 

68  CONTINUE 

69  JB  = JA+1 

Cl  * (GG<  1 , J)-ETA(  JA))/(ETA(  JB)-ETA(JA)) 

DO  70  JJJ=1,4 

70  F(JJJ,J)  = GUJJ»JA)  + (G(JJJ»JB)-G(JJJ,JA))*C1 
DO  71  J=JK, JU 

ETA(J)  = GG(1,J) 

ET(J)  = .5*<ETA(  J)  + ETA(J-1>  ) 

GG ( 3 , J ) = GE/(1.+GE*ET(  J)F 
DO  71  JJ  = 1 , 4 

71  G( JJ,J)  * F ( JJ, J) 

JOEL  = JEO 

G( 2, JDEL)  = 1. 

JK  * JDEL+1 
DO  72  J=JK, JU 
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72 


73 


G(  1 , J-i)*HS(  J)  *G  (2  , J-l > 

1. 

0. 

1 . - UE  P * *2 


.GE.  G ( 2 
JMA  X = J 
.LE.  G ( 2 


J-l)  .AND.  G(2,J)  .GE.  G(2,J+1)  .AND.  JMAX 
J-l)  .AND.  G(2,J)  .LE.  G(2,J+1>>  JM  IN=J 


0)  GO  TO  75 


G( 1, J)  = 

G(  2,J)  = 

G(3,J)  = 

G ( 4 » J ) s 
CONTINUE 
JMAX  = JU 
DO  73  J = 2 » JDE  L 
IF  ( G( 2 » J ) “ 

S .EQ.  JU) 

IF  ( G( 2 • J ) 

CONTINUE 
IF  ( NPRNT  .EQ. 

DO  74  J=1,JU 

PRINT  106, ETA  I J)  ,G!1,J)  ,G(2,  J)  ,G  (3  , J ) ,G  <4,  J ) 

CONTINUE 
CONTINUE 

FORMAT  (/.10H  JET  POINT  »I4»6X,4HXI  = , El 3 .6 , 3X , 4HUE  =,E14.7,2X, 
S5HCPW  =,E13.6,2X,5HCPE  =,E13 .6  ,2X  ,10HCP  t POT  ) =,E11.4) 

101  FORMAT  I20X,6HTAUW  = ,E12. 5 ,2X  ,7HDELX  I = , E12. 5,  IX  , 7HDELTA  =,E11.4, 
$2X,6HGMAX  = ,E1 2 . 5 ,2X ,6HGMI  N =,E12.5) 

102  FORMAT  ( 2 OX  ,9HI  TERATI  ON  ,1 4 ,7X  , 8HET  AMAX  =,F9.4,3X, 

$ 8HETAMIN  = ,F9.  4 ,3  X .5HDHW  = , E13 .6 , 2X  ,6HDCPW  =,E14.7) 

103  FORMAT  ( 2 OX  , 8H  JOELTA  =,I4,8X,4HHW  =,E14.7) 

104  FORMAT  ( / ,8X,3HETA  ,13  X , 1HF  ,1  5X  , 1HG , 1 5X , l HH,  14X,  2HCP,  14X,  4HEPS*, 
t 13X,3HTAU) 

105  FORMAT  (7E16.7) 

106  FORMAT  (5E16.7) 

110  FORMAT  (6E13.6) 

114  FORMAT  ( //,  40H  SEPARATED  PROFILE  CALCULATED  AT  ST  AT  ION,  14,  5X, 

$1 2H0N  I TERATION  ,1  4 ,5  X ,5HC  PW  =,E16.7) 

115  FORMAT  ( / / , 42H  GRID  ADJUSTED  FOR  DOWNSTREAM  CALCULATIONS,//) 

RE  TURN 

END 


74 

75 
86 

100 
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1)  JE  = JMI  N 


AND.  IT  .EQ.  1 )C1=1.05*CI*S  (I,N-I)/E3 


SUBROUTINE  EDDY  ( E 1 ,E2 , E3 ,E4 ,E5 , E6, E7 , JMAX, J MIN, J DEL ,LU, L , IT ) 
DIMENSION  S(  4*501 

C0MM0N/HBDBJE/ETA(121>  ,EPS ( 2 ,12 1 ) , G( 4 , 12 1 ) , GG(7 , 121 ) 

N = L 
JE  = JOEL 
IF  ( JMI  N .GT. 

S(  1,N)  = E3 
Cl  = E 1 *E  2 *E  3 
IF  (N  .GT.  1 . 

C 2 = Cl/15. 

C 4 = 0. 

DO  4 J=2 , JMAX 

4 C 4 = C4+.5*(ETA(  J)-ETA(  J-l))  *(G  (2  , J) +G(2  , J-l ) )*  ( 1.- .5*  ( G(  2,  J ) 
$ ♦G(2,J-D)  /G  ( 2 , JMAX)  > 

S<  2»N)  = C 1 *C4 

IF  (N  .GT.  1 .AND.  S(2,N)  .LT.  S(2,N-1))  S ( 2,N)  = S(  2,  N-l ) 

RET  = S(  2 ,N) 

IF  (RET  .LT.  300.)  RET=300. 

C4  = (RET*.  001 ) **2 

CK1  = .44-.  19/(  l.+.49*C4) 

CK3  = 26.+14./II.+C4) 

C 4 = RET/425. -1. 

IF  ( C 4 .LT.  0.)  C4=0. 

C4  = l.-EXP(-.243*C4**.5-.298*C4) 

CK2  = .0168*1.  55/(1. «-.55*C4) 

C 3 = CK1**2*C1 
JWJ  = 2 

IF  ( JMI N .EQ.  1)  JW  J=1 

IF  ( JMI  N .EQ.  1 .AND.  JMAX  .EQ.  JOEL)  JWJ=0 
IF  ( JWJ  .EQ.  0)  GO  TO  8 
C4  = . 5 *(  G ( 2»JMAX)+G(2  ,JE)  ) 

DO  1 J= JMAX  , JOE  L 
IF  (G(2,J)  .LE.  C4)  GO  TO  2 

1 JL  = J+l 

2 SF  = E TA(  JL)-ETA(  JMAX) 

GO  * G(  2 , JMAX)-G  ( 2 • JF  ) 

GL  = G(  2,  JE)  /GD 

IF  TgL  .LT.  3.04  .ANO.  GL  .GT.  1.1)  SL=.403  + .035*  (GL-1.1 ) 

IF  (GL  .LE.  l.i)  SL*. 541 *(GL-.52) *( G L-l . 1 ) - 1 .1 14*GL* ( GL- 1.1 ) 

S ♦ •632*GL*(GL-.52J 
SL  * SL*SF 
EB  = 1.4 

IF  (GL  .LT.  3.04  .AND.  GL  .GT.  1.1)  EB= l .59-.098* (GL- 1 . 1 ) 

IF  (GL  .LE.  1.11  EB=3.65*(GL-.52)*(GL-i.l)-5.71*GL*(GL-l.l) 
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t +2.49*GL*(GL-.52) 

EB  = ETA(  JMAX)+E9*SF 
S(  3 »N ) = C2*GD*SL 
C4  = E 1 *E  TA  ( JL)  /E6 
CC1  = l.34*C4**2 

CC  = ( 1. +4. 48*C4-CC1) **2/(1. +0011/(1.  +.3*C4) 

S(  3.N)  = S(  3,N)  *CC 

IF  (N  .GT.  1)  S(3,N)=.5*(S(3  ,N)+S(3t  N-l)  ) 

S14.N)  = S(  3 »N) 

IF  (N  .GT.  I)  GO  TO  5 
DO  3 J=2,LU 
JJ  = LU+l-J 

IF  (EPS(lfJJ)  .LT.  EPSUiJJ+lM  GO  TO  5 
3 EPEB  = EPS(  lfJJ+ll-l. 

5 CA  = l.  + S(  4,N) 

C 5 = .5*S(4,N) 

IF  (JWJ  .EQ.  1)  GO  TO  9 
C6  = E TA(  JMAX ) *2.3 *SF 
DO  6 J=JL  »LU 

IF  (ETA(J)  .GT.  C6)  GO  TO  7 

6 JM  = J 

7 C 5 = . 5*E  PEB 

EB  = . 8*E  TA ( JOEL) 

SL  = . 127*E  TA(  JOEL  ) 

GO  TO  9 

8 C 5 = . 5*CK2*C1*(ETA(  JDEL)  *G  ( 2 , JDEL) -G<1  , JDELM 
EB  = . 8*E  TA  ( JDE  L ) 

$L  * . 1 2 7 *E  T A ( JDEL) 

C4  * l.  + 2.*C5 

9 SL  = 1.A1A+SL 

EPO  = l.+CK2*Cl*(ETA( JMAX) *G (2.JMAX )-G( 1 ,JMAX ) ) 

IF  (N  .GT.  I .AND.  IT  .EQ.  1)  Cl=.95  *E1  *E2*E3 
C 6 = (Cl  *G(  3 » I ) I **.5 
C 7 = 1I.8/C6 

C8  = l.+El**2*E2*E4*C7*(l.+.5*El*C7/E6)/  (E3*G(3»  1)) 
IF  ( C 8 .LT.  .7)  C 8=.  7 
C 8 = C8**.5*C6/(CK3*(1.+E1*C7/E6) ) 

DO  10  J=l,JMAX 
C 6 = C 8*E  TA  ( J ) 

C 6 = CKl*ETA( J) *(1.-EXP(-C6)I 
EPS(1»J)  = 1.+C1  *G(3  ,J)  *C6**2 
IF  (EPS(1»J)  .GT.  EPO)  EPSC I » J)  =E  PO 
IF  (EPS(ltJ)  .GE.  EPO)  GO  TO  11 

10  JC  = J 

11  CONTINUE 
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JC  = JCU 

C 6 = .5*(C4-EPS(  1,  JC)) 

Cl  - • 1 75*E  TA(  JMA  X) 

C 8 * . 65*E  TA( JMA  X) 

00  12  J-JC  « JMA X 

12  EPS(l.J)  = EPS(  1 » JC)+C6*(  l.*ERF  ( (ETA<  J) -C8  )/C7  ) ) 
IF  (JWJ  .EQ.  i»  JM= JMAX 

13  IF  (JWJ  .EQ.  0)  JM= JC 
00  14  J=JM»  LU 

14  EPSdtJl  = l.+C5*(l.-ERF( <ETA( J)-EB)/SL) ) 

IF  (JWJ  .LE.  1)  GO  TO  17 

00  15  J = JMAX»JL 

15  EP  S(  1*  J)  = C4 

C 5 = EPS(  1,JM)-C4 
C 6 = 1.570796/<ETA(  JM)-ETA(  JL)  ) 


DO  16  J=>JL  » JM 

16  EPS(1,J)  * C4+C5  *( 1.  -COS (C6*  ( ET  A ( J ) -ET  A( JL ) ) ) ) 

17  DO  18  J=1  t LU 

18  EP  S(  2 » J ) = EPS(  l ,J)  *GG(4,  J) 

C 7 = -.1387/SF 

C 6 = E TA(  JMA X)  85*SF 
C 5 = .6 

C 4 = l.  /C5**3/ETA(  JMAX) 

C5  = C 5*E  TA ( JMAX ) 

DO  20  J=1 »LU 
C 8 = 0. 

IF  (ETA(J)  .LT.  C 5)  C8=C4*ETA(J) 

IF  (ETA(J)  .GE.  C5  .AND.  J . LE.  JMAX)  C8=  ( ET  A ( JM  AX  )/  ET  A(  J ) )**  2 
IF  (J  .GT.  JMAX  .AND.  ETA(J)  . LE.  C6 ) C8  = l. -1.387* 

S ( (ETA(  J)-ETA(  JMAX)  )/SF)  **2 
20  GG(6,J)  = C7*C8 
RE  TURN 
END 
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SUBROUTINE  RADIUS  <X»RC) 
KIND  AIRFOIL 
RC  = .037777778 
IF  ( X .LT.  . 001  I RC  = • 1 
IF  ( X .EQ.  0.  ) RC  =•  408 
RE  TURN 
END 


7 


f 
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SUBROUTINE  INVERT  (D) 

DIMENSION  D(4,4)  ,U(4  ,4)  ,UI(4,4) 

INVERT  A 4X4  MATRIX 
DO  1 J=I,4 

1 U( J,I)  = D ( J,l) 

DO  2 J=2,4 

U< i,J)  = DC  l ,J)  /U(  1 ,1) 

2 U ( J , 2 ) = DC  J,2)-U(  J,i)  *U(  1 ,2) 

DO  3 J=3  r4 

U(2,J)  = (D(2»J)-U(2*l)  *U(  1 , J)  ) /U (2  ,2 ) 

3 U(  J » 3 ) = DC  J,3)-U(  J,l)*U(i,3)-U(  J,2)*U(2,3) 

U( 3»  4)  = (D(3,4)-U(3,l)*U(l,4)-U(3,2)*U(2,4) )/U<3.3) 

U(4,4)  = D(  4f4)-U(4  ,1 ) *U(1  ,4)  -U(4  ,2)  *U  (2  .4  )-U  (4,  3)*U  (3,  4) 

UKIil)  = UC2.2)  *U(3,3)  *U(4,4> 

UI(2,1>  = -U(2,i)*U(3.3)*U(4,4> 

UI  ( 3 * 1 ) = U(2,1)*U(3,2)*U(4,4)-U(3,1)*U(2,2)*U<4,4> 

UI(4,1)  = U(2fl)  *U(4»2)  *U(3,3)*-U(3»i)*U(2,2)*U(4,3) 
i -U(2,1)*U(  3,2)*U(4,3)-U(2,2)*U(3,3)*U(4,  l) 

UI  ( 1,2)  = -U(l  .2) 

UI  ( 2,2)  = U(I.I)  *U(3,3)  *U(4,4> 

UI  ( 3,2)  = -U(l,l)*U(3»2)  *U(4  ,4) 

UI  C 4 » 2)  = U(  l,l)*U(3,2)*U(4,3)-U(l,l)*U(4,2)*U(3,3) 

UI C 1,3)  = U ( I v 2 ) *U  ( 2 ,3 ) -U  ( 1 , 3 ) 

UI(  2,3)  = -U(2,3) 

UI  ( 3 » 3 ) = UI  1 « 1 1 *U(2  ,2)  *U(4,4) 

UI(  4 » 3 ) = -U<  I ,1)*U<4,3)  *U(2  ,2) 

UI(  1,4)  = -U(1,2)*U12,3)  *U(3  ,4 ) +U  ( l , 2 )*U  ( 2, 4 ) +U  ( 1,  3 )*U  ( 3,  4)-U(  1,4) 
UK  2 ,4)  = U(  2 ,3)  *U(  3 ,4)  -U(2  ,4) 

UK  3,4)  = - U(  3 ,4) 

UI( 4 ,4)  = U(  1,1) *U( 2 *2) *U(3 ,3) 

DO  4 J=l,4 

4 D(  4,  J)  = UI  (4,  J) 

DO  5 J=1 ,3 

D ( J , 4)  = UI  ( J ,4 ) *UI  (4,4) 

D(  3,  J ) * UI  ( 3,  J)  + UI  (3  ,4)  *UI  (4  , J) 

D ( 2,  J ) = UI(2,3)*UI(3,J)+UI(2,4I*UIC4,J) 

5 D(  1 , J ) = UI(  1,3)  *UI  ( 3.JH-UI  (I,4)*UI  (4,J) 

D(  2,2)  = DC  2,2)  + UI  (2  ,2) 

D(  2,1)  = D(2,1)  + UI(2  ,1) 

DC  1,2)  * DC  1 ,2)+UI(  1 ,2)  *UI  (2  ,2) 

DC  1,1)  * DC  1 tlH-UI  ( 1 ,1M-UI  Cl  ,2)  *UI(2,l) 

C = U(  l,l)*U(2,2)*U(3,3)*U(4,4) 

DO  6 J=1 , 4 
DO  6 K=1 ,4 

6 DC  J»K)  = DC  J,K)/C 
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C = ABS(C) 

IF  (C  .LT.  l.E-3)  PRINT  7.C 
7 FORMAT  ( 38H  NEARLY  SINGULAR  MATRIX, 
RE  TURN 
END 


DETERMINANT 


, E18.7) 
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r - 4 rl.r  • , C 2 = UT/SI3) 

’f^^HSJET  » .0005 

v f & jf  si'up  ^?v»  ji  RS  = 1.13 

JU  = 121 


SUBROUTINE  SLOTS  L 
DIMENSION  A<  5,241)  .ETAI121) 

COMMON /JE  TMI  X/HS(12l)  ,F(6,121),S<18) 
COMMON  /BL  SLB  L/ED  V(  5,241)  » NBL » NEOGE 
COMMON  /MAI  NSL/STFT(/),LAMSL 
EQUIVALENCE  ( A ( 1 ,1 ) ,EDV(  1 ,1 ) ) 

DO  1 J = l,7 
S(  J)  = STF  T ( J) 

CS  = S(  1) 

UT  = S(  6 ) 

Cl  = CS*S(2>  *UT 


HS(  1)  = HSJET/RS 
ETA(l)  = 0. 

DO  2 J=2  , JU 

HS(  J)  = HS(  J-l ) *RS 

IF  ( H S(  J ) .GE.  .04)  HS(  J)  =.04 

ETA(  J ) = ETA(  J-l  H-HSI  J) 

IF  (ETA(J)  .GT.  .006)  RS=l.l 
IF  (ETA(J)  .GT.  .1)  RS=1 • 1 5 
JCL  = J 

IF  (ETA(J)  .GT.  .45)  GO  TO  3 

2 CONTINUE 

3 HCL  = .5-ETAIJCL) 

H S(  JCL*  1 ) = 2.  *HCL 

ETA(  JCL*1 ) = E TA  ( JCL) -mS  ( JCL+1  ) 

DO  4 J=2  , X L 
K = JCL+J 
I = JCL+2-J 
HS(  K)  = H S(  I ) 

4 ETA(K)  = ETA(  K-l  )+HS(K) 

JD  = 2*JCL 

Cl  = 0. 

F ( 1 » 1 ) = 0. 

F(  2,1)  = 0. 

F(  3,1)  = Ci*C2 
F( 5,1)  = l. 

00  5 J=2,  JCL 

F ( 5 , J ) = l.*.4*(Cl*ETA(J)  ) **3/  ( ( Cl *ET  Af  J) ) **2*324. ) 
$ /(  l.*6.  8*C  S*E  TA(  J)  / S(  1 ) ) 

IF  (LAMSL  .EQ.  1)  F(5,J)  =1.  . 

F(  3,J)  = F(  3,1)/F(5,J)*(1.-2.*CS*ETA(J)/S(1)  ) 


179 


F(2,J)  = F(2»J-l)*.5*HS(J)*(F(3,J)  + F(3,J-l)> 
F ( 1 f J ) = F(  1 .J-l  ) + .5*HS(  J)*|FI2,J)»f(2.J-in 

5 Cl  = Cl  ♦.  5*H  5(  J)*(F(2»J)**2*-F(2»  J-l ) **Z  ) 

Cl  = 2.  *C  I + 2.  *HC  L*F  ( 2 » JCL)  **2 

F ( 1 » JCL  + 1 ) = F(  1 ,JCL)+2.  *HCL*F  (2,  JCL) 

DO  7 J=1 » X L 

K = JCL+J 

I = JCL+l-J 

IF  (J  .EQ.  I)  GO  TO  6 

F ( 1 1 K ) = FC 1 ,K-1 ) + FC l,I+l)-F(l,I  ) 

6 F ( 2»K)  = F C 2 * I I 

F ( 3»K)  = -F ( 3 » I ) 

7 F ( 5»K)  = F ( 5 » I ) 

CMU  * 2.  *C  S*S(  3)  **2*CI 
UAVE  = C S*S(  3)  *F  ( l » JD) /S ( 1 ) 

CJ  = CS*S(  3 ) *F  ( 1 t JD ) 

UMAX  = S(  3)  *F(  2 , JCL) 

PRINT  IOI 

PRINT  I02»CMU,C  J ,UAVE tUMAX »F (3  ,1) 

JBL  = NEDGE+10 

IF  (JBL  .GT.  NBL ) JBL  = NBL 

JLO  * JU-JD 

JO  = ( 6*JBL-3*JL0) /(2*JLO)*-l 
JR  = JLO/3 
JL  = J3  + 1 
JR  1 = JR*- JL 
JR  2 = 2*JR+ JL 

C 3 = ( 2. *S(  5))**.5/(CS*S(2)**.  5*S(3)  ) 

K = 1 
LL  = 0 

DO  8 J=JL,JU 
LL  = LL+1 

IF  (LL  .EQ.  I .AND.  J . LE . JRil  GO  TO  8 
K = M-l 

IF  (J  .GT.  JR2 ) K*K+J0 
ETA(J)  = ETA(  JD)+C3*A(1  ,K) 

F(1,J)  = F(  1,  JD)+C3*A(2,K) 

F ( 2,J)  = A ( 3 » K) 

F ( 3,J)  = A(  4,K)  /C3 
F ( 5,J)  = A( 5.K) 

8 IF  (LL  .EQ.  I)  LL=-1 
DO  9 J=JL  » JR1 *2 

ETA(J)  = .5*(ETA(  J-1)+ETA(  J*l)  ) 

F ( 1 * J ) = .5*(F( L,J-l)+F(l»J+l>) 

F(2,J)  = .5*(F(2tJ-I)+F(2tJ+I)) 
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F ( 3 » J ) = • 5 *(  F ( 3 » J- 1 ) +F  ( 3 t J*  1 ) ) 

9 F ( 5 » J ) = .5*(F(5,J-1)  + F(5  , J+l)  ) 

00  11  J=1,JU 

IF  (J  .GT.  JO)  HSI  J)  =ETA(  J) -ETA(  J-ll 

11  CONTINUE 
S(  5)  = 0. 

DO  12  J=1,JU 

IF  (J  .GT.  JL  .AND.  F(2,J)  .GT.  .99  .AND.  S(5)  .EQ.  0.) 

S S(  5)  = S(  1 ) *E  TA  ( J) 

12  F(4.J)  * S<4) 

S(  4 1 = S(  3) 

NPRNT  = 0 

IF  (NPRNT  .EQ.  0)  GO  TO  13 
DO  10  J = 1 , JU 

10  PRINT  104,ETA(  J)  ,F(1  ,J)  ,F(2,  J)  ,F(3,J»  ,F(4,  J ) , F(  5,  J ) 

13  CONTINUE 

101  FORMAT  ( /// » 10X  »24H  JET  STARTING  CONDITIONS) 

102  FORMAT  {/,6H  CMU  =fF8.  5 ,5X  ,4HCM  =,F  9 .6 ,5X  ,6HU  AV  E =,F7.3»5Xf 
S 6HUMAX  =,F7.3,5X,4HHW  =,F7.2) 

104  FORMAT  (6E13.6) 

RE  TURN 
END 
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